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Abstract 

A  numerically  stable,  accurate,  and  robust  form  of  the  exponential 
characteristic  (EC)  method,  used  to  solve  the  time-independent  linearized 
Boltzmann  Transport  Equation,  is  derived  using  direct  affine  coordinate 
transformations  on  unstructured  meshes  of  tetrahedra.  This  quadrature,  as 
well  as  the  linear  characteristic  (LC)  spatial  quadrature,  is  implemented  in 
our  transport  code,  called  TETRAN.  This  code  solves  multi-group  neutral 
particle  transport  problems  with  anisotropic  scattering  and  was  parallelized 
using  High  Performance  Fortran  and  angular  domain  decomposition.  A  new, 
parallel  algorithm  for  updating  the  scattering  source  is  introduced.  The  EC 
source  and  inflow  flux  coefficients  are  efficiently  evaluated  using  Broyden's 
rootsolver,  started  with  special  approximations  developed  here.  TETRAN 
showed  robustness,  stability  and  accuracy  on  a  variety  of  challenging  test 
problems.  Parallel  speed-up  was  observed  as  the  number  of  processors  was 
increased  using  an  IBM  SP  computer  system. 
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DEVELOPMENT  OF  A  DISCRETE  ORDINATES 
CODE  SYSTEM  FOR  UNSTRUCTURED  MESHES 
OF  TETRAHEDRAL  CELLS,  WITH  SERIAL 
AND  PARALLEL  IMPLEMENTATIONS 

Chapter  I:  Introduction 

We  present  the  development  of  a  parallel  unstructured  tetrahedra 
mesh  discrete  ordinates  radiation  transport  code.  This  code,  TETRAN,  solves 
the  linear  time-independent  Boltzmann  transport  equation  for  the  angular 
flux  density  using  the  linear  and  exponential  characteristic  spatial 
quadratures.  It  also  uses  the  standard  source  iteration  algorithm  with  multi¬ 
group  cross  sections  and  general  anisotropic  scattering.  The  code  was  made 
parallel  using  High  Performance  Fortran  directives.  We  examine  the 
performance  of  the  code  using  several  benchmark  problems  and  present  the 
results.  The  exponential  characteristic  method  was  found  to  be  a  stable, 
robust  coarse  mesh  method  when  provided  with  positive  source  moments. 
Several  challenging  test  problems  demonstrate  EC's  performance. 
Performance  of  the  parallel  code  is  briefly  presented. 

Background 

The  fundamental  equation  describing  neutral  particle  transport 
processes  is  the  linear  Boltzmann  transport  equation  (BTE)  (Lewis  and 
Miller,  1993:1).  The  time-independent  BTE  describes  the  balance  of  neutral 
particles  in  a  six-dimensional  phase  space:  three  spatial  dimensions  (r),  two 


angular  dimensions  (Q),  and  energy  (E).  The  BTE  relates  the  rate  of  change 
of  neutral  particle  density  at  a  point  r ,  moving  with  energy  E,  in  direction  Q , 
to  the  differences  between  sources  and  losses: 

[Q  •  V  +  a(r  ,E)]v|/(r,Q,E)  =  S(r,Q,E).  (1) 

where  the  angular  flux  density,  \\i ,  is  the  flux  of  particles  at  point  r,  moving 
with  energy  E ,  in  direction  Q .  The  bracketed  term  accounts  for  particle 
losses  due  to  streaming  and  collisions.  S(r,Q,E)  represents  the  source  of  all 
particles  in  the  (r,Q,E)  phase  space  cell.  The  characteristic  solution  to  the 
BTE  is  presented  in  Chapter  2.  In  Chapter  3,  we  discuss  the  multi-group 
energy  discretization  of  the  BTE  and  the  impact  of  the  traditional  anisotropic 
scatter  treatment  on  the  parallelization  of  the  source  iteration  algorithm. 

Analytical  solutions  for  the  BTE  exist  for  idealized  problems  (one 
energy  group  with  isotropic  scatter  and  symmetry  boundaries).  Recently, 
researchers  at  Los  Alamos  developed  several  semi-analytic  benchmarks  that 
solve  multi-dimensional  transport  problems  (Kornreich,  1997).  Point  kernel 
techniques  are  used  in  some  radiation  shielding  calculations  where  speed, 
rather  than  accuracy,  is  required  to  determine  dose  rates  from  simple  source 
configurations  (Glasstone,  1981).  Build-up  factors,  calculated  via  the 
numerical  techniques  discussed  below,  often  augment  point  kernel 
calculations. 
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Numerical  methods  are  generally  used  to  solve  the  BTE.  The  two  most 
popular  classes  of  methods  are  Monte  Carlo  and  deterministic.  Where 
diffusion  processes  dominate  the  transport  of  particles,  simplifying 
assumptions  can  be  made  to  the  BTE,  transforming  it  into  the  diffusion 
equation  which  is  solved  by  well  established  algorithms  (such  as  finite 
difference  or  conjugate  gradient)  (Ott,  1983). 

Monte  Carlo  Techniques 

The  Monte  Carlo  technique  simulates  individual  particle  path  histories 
and  accumulates  the  results.  As  these  histories  are  accumulated,  statistics 
are  formed  that  estimate  the  uncertainties  in  the  results.  The  method  is 
generally  suitable  to  complex  three-dimensional  geometries.  Computational 
grids  are  not  explicitly  needed.  Besides  its  general  geometry  capability,  the 
Monte  Carlo  technique  can  also  take  advantage  of  continuous  energy  nuclear 
data  which  is  the  natural  form  of  cross  section  data  found  in  the  Evaluated 
Nuclear  Data  Files  (ENDF)  (Rose,  1990).  It  can  also  use  multigroup  cross 
section  data  that  are  typically  used  for  deterministic  transport  algorithms 
(Briesmeister,  1991).  The  primary  disadvantage  of  the  technique  is  its  need 
for  variance  reduction  techniques  for  complicated,  deep-penetration 
problems.  Another  disadvantage  is  that  the  method  can  not  represent  the 
flux  in  an  entire  system  except  at  great  expense  since  numerous  flux 
estimators  are  required.  Deterministic  methods  directly  create  a  flux  map. 

An  excellent  discussion  of  the  Monte  Carlo  method  is  found  in  Lewis  and 
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Miller,  Chapter  7  (Lewis  and  Miller,  1993),  or  the  MCNP™  4B  (Monte  Carlo 
N-Particle  Transport  Code  System)  manual,  Chapter  1  (Briesmeister,  1991). 
MCNP™  is  widely  used  due  to  its  general  three-dimensional  geometry 
ability,  accuracy,  and  ease  of  use;  we  use  it  here  to  generate  benchmarks  for 
comparison  with  TETRAN  results. 

Deterministic  Methods 

Deterministic  methods  solve  the  BTE  using  simplifying  assumptions 
and  techniques  to  perform  the  angular  and  spatial  integrations.  The  most 
widely  used  method  to  solve  the  BTE  is  discrete  ordinates  (Lewis  and 
Miller,  1993:  116).  First  utilized  for  transport  calculations  in  stellar 
atmospheres  (Chandrasekhar,  1960),  researchers  at  Los  Alamos  National 
Laboratory  adapted  this  method  for  neutron  transport  work  in  the  late  1950's 
(Carlson:  1958). 

In  the  discrete  ordinates  method,  the  BTE  is  solved  for  a  finite  number 
of  angular  directions,  Qn ,  to  get  the  angular  flux,  v|/n  ,  for  each  cell  in  a 
computational  mesh.  After  all  of  the  vj/  n  have  been  computed  for  a  given 

mesh,  numerical  quadrature  is  used  to  calculate  the  meaningful  integral 
quantity:  the  scalar  flux,  <j) ,  in  each  cell  i, 

N 

4>i  =  ZWnVi,n-  (2) 

n=l 
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Usually,  the  angles,  Qn  ,  and  weights,  wn ,  are  chosen  to  accurately  integrate 
spherical  harmonics  functions,  as  is  described  in  detail  in  Lewis  and  Miller, 
Chapters  3  and  4  (Lewis  and  Miller,  1993). 

Many  algorithms  exist  within  the  discrete  ordinates  family  of  methods. 
One  of  the  oldest,  yet  still  used,  is  the  diamond  difference  method  (Lewis  and 
Miller,  128-131).  The  diamond  difference  method  makes  a  finite  difference 
approximation  to  the  differential  operator  in  the  BTE  and  assumes  several 
auxiliary  equations.  The  method  is  easy  to  implement  and  is  computationally 
inexpensive.  However,  diamond  difference  suffers  from  numerical 
instabilities  that  produce  oscillations  in  the  angular  flux  (Petrovic,  1996)  and 
generally  does  not  propagate  flux  in  the  correct  direction  (Mathews,  1998).  In 
slab  geometry,  diamond  difference  can  also  produce  negative  angular  fluxes 
unless  the  computational  cell  optical  thickness,  s ,  is  small  (<  0.1).  In  two- 
and  three-dimensions,  diamond  difference  is  not  even  conditionally  positive 
unless  fix-ups  are  used.  This  behavior  limits  the  application  of  diamond 
difference  to  deep-penetration  shielding  problems  since  exceedingly  fine 
meshes  are  required  to  maintain  the  positivity  and  accuracy  of  the  method. 
Indeed,  the  diamond  difference  method  or  one  of  its  cousins  [adaptive 
weighted  diamond  difference  (AWDD)  (Alcouffe,  1993)  and  theta-weighted 
diamond  difference  (Rhoades  and  Azmy,  1996)]  has  been  the  backbone  of 
every  major  discrete  ordinates  neutral  particle  transport  code  developed  in 
the  United  States  since  the  inception  of  the  method  around  1960. 
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Figure  1  shows  the  progression  of  production  codes,  the  methods  they 
employ,  and  a  backdrop  of  computational  hardware  generally  available  at  the 
time  of  the  method  [(Intel  Corp.,  1998)  (Foster,  1995),  (Dongarra,  1998), 
(USDOE  ASCI,  1998),  (Wareing,  1997),  (Lewis  and  Miller,  1993),  (Rhoades, 
1995),  (Alcouffe,  1995),  (Briesmeister,  1991)].  The  figure  shows  that  few  new 
(higher  order)  spatial  quadratures  have  been  incorporated  into  the 
production  codes  in  the  last  30  years.  The  reason  for  this  is  that  higher  order 
methods  are  generally  more  computationally  intensive  per  spatial  cell  than 
the  traditional  or  weighted  difference  approaches.  This  is  particularly  true  of 
the  characteristic  methods  presented  in  this  dissertation.  Much  of  the 
apprehension  in  embracing  higher  order  methods  stems  from  Lathrop's 
seminal  paper,  "Spatial  Differencing  of  the  Transport  Equation:  Positivity  vs. 
Accuracy"  (Lathrop,  1969).  In  this  paper,  Lathrop  presents  what  have  become 
the  de  facto  criteria  for  a  "desirable"  spatial  quadrature: 

1.  "it  should  be  accurate  in  the  sense  that  it  has  small  truncation  error, 

2.  it  should  be  simple,  which  means  that  it  should  involve  a  small  number  of 
numerical  operations  and  should  involve  unknowns  within  a  single  mesh 
cell, 

3.  it  should  produce  positive  fluxes  if  the  source  and  boundary  fluxes  are 
positive, 

4.  the  scheme  should  be  conservative  in  the  sense  that  a  well  defined 
relation  exists  among  flows  into  and  out  of  the  cell  and  sources  and 
absorptions  in  the  cell,  and, 

5.  the  scheme  is  easily  generalized  to  all  (curvilinear)  geometries." 
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Obviously,  the  above  criteria  are  reasonable,  particularly  from  the 
vantage  point  of  1969  when  computational  hardware  constraints  were  severe. 
However,  the  2nd  and  5th  criteria  have  prevented  the  acceptance  of  many 
higher  order  methods  because  they  tend  to  be  computationally  expensive  per 
cell,  not  very  simple,  and  are  not  generalizable  to  curvilinear  geometries. 
Criterion  5  proves  to  be  a  particular  problem  for  characteristic  methods  since 
such  methods  assume  particle  flow  along  straight  lines  in  Cartesian  space. 
However,  we  argue  that  criteria  2  and  5  are  unreasonable  given  the 
computational  and  mesh  generating  capabilities  that  exist  today.  Indeed,  for 
true  three-dimensional  problems  not  involving  the  types  of  symmetries 
needed  to  exploit  curvilinear  coordinates,  an  unstructured  Cartesian  mesh  is 
probably  a  better  choice  for  performing  transport  analysis. 
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Correlation  of  Production  Code  Transport  Algorithms 
to  Hardware  Floating-Point  Performance 
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Figure  1.  Correlation  of  Production  Code  Transport  Algorithms  to  Hardware 

Floating-Point  Performance. 

Figure  1  shows  the  fantastic  growth  in  floating-point  performance  from  the 
days  of  the  first  computer  to  near  term  plans  for  the  next  generation  of 
supercomputers  envisaged  by  the  Department  of  Energy's  Advanced 
Strategic  Computing  Initiative  (ASCI).  Clearly,  system  performance  has 
grown  tremendously  since  the  birth  of  computational  radiation  transport 
circa  1960  (both  in  processor  performance  and  memory).  Meanwhile, 
transport  researchers  have  derived  more  accurate  transport  methods  than 
the  traditional  or  weighted  diamond  difference  approaches  which  are  still  the 
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primary  methods  of  choice  for  the  large  discrete  ordinates  production  codes: 
DANTSYS  (Alcouffe,  1995)  and  DOORS  (RSICC,  1998).  Additionally,  most  of 
the  modern  discrete  ordinates  transport  codes  use  Cartesian  computational 
meshes  with  boxoid  (in  the  case  of  three-dimensional  transport)  or 
rectangular  (in  the  case  of  two-dimensional  transport)  computational  cells. 
This  approach  simplifies  meshing  a  problem  where  curved  surfaces  are  not 
involved  or  one  does  not  want  to  resolve  fine  details  of  the  problem.  However, 
the  number  of  cells  needed  to  mesh  a  problem  accurately  grows  quickly  with 
the  detail  desired  to  be  refined.  It  is  for  this  reason  that  researchers  continue 
to  pursue  the  development  of  higher  order,  linear  and  non-linear  discrete 
ordinates  methods  for  use  on  unstructured  grids. 

Non-linear  Discrete  Ordinates  Methods 

Over  the  past  ten  years,  researchers  at  the  Air  Force  Institute  of 
Technology  (AFIT),  Texas  A&M  University,  and  Los  Alamos  National 
Laboratory  (LANL)  have  developed  several  non-linear  characteristic  methods 
to  solve  the  BTE.  In  these  non-linear  methods,  the  needed  source  or  flux 
coefficients  are  obtained  by  root-solving  a  non-linear  system  of  equations. 
These  methods  were  pursued  due  to  their  inherent  positivity  given  positive 
source  and  boundary  data  and  their  potential  for  robust  coarse  mesh 
performance.  The  timeline  showing  the  development  of  non-linear  methods 
research  is  shown  in  Figure  2. 
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NLC:  Non-Linear  Characteristic 
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Figure  2.  Development  timeline  of  non-linear  methods  research. 

First  among  this  family  of  methods  were  the  step  adaptive  (Mathews, 
1990)  and  linear  adaptive  (Mathews,  1993)  methods.  Following  the 
development  of  the  adaptive  methods  was  the  exponential  characteristic 
(EC)  method  for  slab  geometry,  developed  at  both  AFIT  (termed  exponential 
characteristic  or  EC)(Mathews,  1994)  and  LANL  (termed  non-linear 
characteristic  or  NLC)(Walters,  1996).  This  method  takes  the  average  and 
first  spatial  moments  of  the  characteristic  solution  to  the  BTE  and  an 
assumed  exponential  distribution  of  the  scattering  source  in  a  cell  to 
calculate  average  and  first  moments  of  the  angular  flux.  Both 
implementations  exhibited  fourth  order  convergence  and  excellent  deep 
penetration  performance  warranting  further  investigation  of  the  method  in 
higher  dimensions.  The  AFIT  implementation  differed  substantially  from 
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LANL's  approach  in  that  the  focus  of  the  development  was  on  producing  a 
numerically  robust  method.  The  LANL  approach  had  obvious  numerical 
conditioning  problems.  Both  groups  continued  the  development  of  EC/NLC  by 
implementing  the  method  in  two-dimensions.  AFIT  developed  both  Cartesian 
rectangular  cell  (Minor,  1995)  and  unstructured  triangular  mesh  (Mathews, 
1997)  EC  algorithms.  LANL  researchers  produced  a  Cartesian  rectangular 
cell  algorithm  (W alters,  1995)  and  convergence  acceleration  techniques  for 
their  methods  in  slab  and  XY  geometry  (Wareing,  1994)  (Wareing,  1995). 
AFIT  continued  our  development  by  extending  the  method  to  three- 
dimensions,  implementing  EC  on  grids  of  unstructured  tetrahedra  (Brennan, 
1996).  The  LANL  researchers  chose  to  pursue  different  methods  (due  to  the 
high  computational  cost  per  phase  space  cell  of  the  NLC  method),  beginning 
the  development  of  the  exponential  discontinuous  (ED)  method  which  they 
have  implemented  in  slab,  two-  (Wareing,  1995)  and  three-dimensional 
Cartesian  cells  (Wareing,  1996).  At  Texas  A&M,  researchers  developed  the 
non-linear  corner  balance  (NLCB)  method  which  they  have  derived  for  slab 
and  two-dimensional  arbitrary  grids  (Castrianni,  1998).  Additionally,  this 
group  also  developed  convergence  acceleration  techniques  (transport 
synthetic  acceleration)  for  use  with  their  method  (Ramone,  1997)  (Anistratov, 
1996). 
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Other  Unstructured  Grid  Transport  Methods 

Besides  the  above  non-linear  methods  for  two-  and  three-dimensional 
transport  on  unstructured  grids,  researchers  developed  several  linear 
methods.  LANL  pioneered  the  development  of  the  linear  characteristic 
method  for  unstructured  triangles  by  developing  a  linear  characteristic-nodal 
method  for  such  grids  (Paternoster,  1984).  A  version  of  the  step  characteristic 
method  was  extended  to  two-dimensional  arbitrary  grids  by  researchers  at 
Oak  Ridge  National  Laboratory  and  the  Westinghouse  Savannah  River 
Laboratory  (DeHart,  1994).  LANL  developed  —  in  concert  with  several  oil 
and  natural  gas  companies  —  a  three-dimensional  unstructured  tetrahedral 
mesh  discrete  ordinates  transport  code.  This  code,  ATTILA,  performs  particle 
transport  using  the  linear  discontinuous  method  with  diffusion  synthetic 
acceleration  of  the  source  iteration  (Wareing,  1996)  and  is  third  order 
convergent  (Wareing,  1998). ).  At  AFIT,  the  step  and  linear  characteristic 
methods  were  developed  for  unstructured  triangular  (Miller,  1993)  and 
Cartesian  rectangular  (Minor,  1993)  cell  meshes.  Brennan  included  a 
development  of  both  the  step  characteristic  and  linear  characteristic  methods 
for  tetrahedral  cells  (Brennan,  1996). 

Purpose 

The  purpose  of  this  research  is  to  contribute  to  the  development  of  a 
high  order,  high  fidelity  (with  respect  to  geometry),  three-dimensional 
discrete  ordinates  transport  code  system.  Currently,  only  one  such  code  exists 
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—  LANL's  ATTILA  (although  it  is  not  a  general  distribution  code).  ORNL's 
TORT  code  uses  the  linear  characteristic  spatial  quadrature  (Rhoades,  1995) 
for  deep  penetration  problems.  However,  it  performs  transport  on  Cartesian 
grids  and  is  thus  poor  at  modeling  geometry  in  detail.  MCNP™,  a  Monte 
Carlo  code,  is  widely  used  for  high  fidelity,  three-dimensional  transport 
calculations.  We  seek  to  develop  a  discrete  ordinates  code  system  that 
performs  neutral  particle  transport  on  arbitrarily  complex  three-dimensional 
geometries  with  accuracy  approaching  or  exceeding  that  of  MCNP™  or  other 
Monte  Carlo  methods.  We  envision  that  such  a  code  system  will  fill  the  gap 
between  Monte  Carlo  and  discrete  ordinates  transport  codes,  producing 
accurate  high  fidelity  global  flux  distributions  for  coarse  meshes  and  deep 
penetration  problems. 

Goals 

The  goal  of  this  research  is  to  implement  the  linear  and  exponential 
characteristic  methods  in  a  parallel  radiation  transport  code.  Specifically,  the 
exponential  characteristic  method  (tetrahedron  computational  cell)  will  be 
derived  using  the  direct  affine  transformation  approach  (Mathews,  1998)  for 
tetrahedral  cell  splitting  and  moment  passing.  This  development  improves  on 
that  of  prior  work  (Brennan,  1996).  Additionally,  we  will  implement  the  EC 
spatial  quadrature  within  a  code  that  is  parallel  capable  with  the 
introduction  of  High  Performance  Fortran  (HPF)  directives. 
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We  begin  this  development  with  the  long-term  vision  of  fielding  a 
production  code  system.  With  that  end  in  mind,  TETRAN  will  be  able  to 
perform  multi-group,  anisotropic  scattering  problems  in  a  high-performance 
parallel-computing  environment.  TETRAN  will  operate  on  arbitrarily 
complicated  meshes  (limited  by  the  mesh  generator  and  available  hardware) 
and  read  standard  cross  section  libraries.  Additionally,  it  will  produce  data  in 
a  format  that  is  easily  read  by  state  of  the  art  scientific  visualization  software 
such  as  AVS/Express. 

Scope 

TETRAN  was  created  to  be  the  kernel  of  a  general  three-dimensional 
neutral  particle  transport  code.  It  solves  the  time-independent  multi-group 
BTE  using  the  EC  and  LC  spatial  quadratures  on  meshes  of  tetrahedra. 
Anisotropic  scattering  is  treated  using  the  traditional  spherical  harmonics 
approach.  The  code  operates  in  parallel  on  an  IBM  SP  when  compiled  with 
the  PGHPF  2.4  development  compiler. 

Limitations 

As  a  first  generation  code,  TETRAN  has  several  limitations  due  its 
maturity  level.  These  limitations  are: 

1)  No  Convergence  Acceleration.  Adapting  or  developing  convergence 
acceleration  may  prove  to  be  straightforward,  but  issues  of  efficiency, 
effectiveness,  robustness,  and  implementation  make  this  suitable  as  the 
focus  of  a  follow-on  effort. 

2)  No  fission  source. 
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3)  Downscatter  only.  To  be  useful,  both  2)  and  3)  would  need  convergence 
acceleration.  Except  for  a  few  unusual  applications,  adding  these  features 
will  be  straightforward. 

4)  High  Performance  Fortran  (HPF).  TETRAN  has  HPF  directives  and  has 
been  minimally  tested  using  them  but  current  compilers  are  still  under 
development. 

5)  Only  vacuum  boundaries  and  no  re-entrant  geometries.  Other  boundaries 
would  unduly  complicate  parallelism  for  a  first  code. 

6)  No  external  incident  beams  or  fluxes.  Same  as  5). 

Several  problems  were  devised  to  test  the  operation  of  the  code  and  it's 
underlying  algorithms.  Each  problem  uses  an  isotropic  emitter  uniformly 
distributed  through  the  volume  of  one  region  as  the  only  source  of  particles  in 
the  problem.  No  incident  current  problems  were  tested  since  the  author  did 
not  know  how  to  use  MCNP™  for  this  type  of  problem.  Additionally,  no 
investigation  of  TETRAN  performance  on  thick,  diffusive  problems  was 
performed,  because  such  problems  require  convergence  acceleration,  which  is 
a  feature  we  have  yet  to  develop  or  adapt  to  these  meshes. 

Organization 

The  remainder  of  this  document  is  composed  of  five  chapters.  Chapter 
2  presents  the  derivation  of  the  exponential  characteristic  spatial  quadrature 
using  direct  affine,  transformations.  Chapter  3  outlines  the  needed  multi¬ 
group  and  anisotropic  scattering  theory  and  presents  the  pseudo-code  for 
TETRAN  as  well  as  the  parallelization  strategy  using  HPF.  We  present  the 
results  of  our  testing  in  Chapter  4.  Finally,  Chapter  5  presents  our 
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conclusions  and  recommendations  for  future  work.  Four  appendices  are 
included.  In  Appendix  A,  we  present  several  exponential  moment  function 
identities  used  in  this  work.  Appendix  B  contains  the  algorithms  used  to  find 
the  source  and  inflow  face  flux  coefficients  needed  for  the  EC  method. 
Appendix  C  summarizes  the  robust  algorithms  used  to  calculate  the  EC 
spatial  quadrature.  Finally,  Appendix  D  presents  the  derivation  of  a  new 
spatial  quadrature,  the  surface  cell  algorithm,  that  was  derived  but  not 
implemented  in  this  effort  because  no  mesh  generators  currently  support  its 
inclusion  in  finite  element  unstructured  meshes. 
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Chapter  II:  Derivation  of  Exponential  Characteristic  Spatial 

Quadrature 

In  this  chapter,  we  derive  the  exponential  characteristic  (EC)  solution 
to  the  Boltzmann  Transport  Equation  (BTE)  for  a  tetrahedron  spatial  cell. 
Specifically,  the  equations  for  the  average  and  first  spatial  moments  of  the 
angular  flux  in  a  tetrahedron  cell  are  derived  as  well  as  the  needed  source 
and  inflow  flux  systems  of  equations.  The  chapter  ends  with  a  brief  summary 
of  the  linear  characteristic  method,  which  was  implemented  prior  to  the  EC 
quadrature  in  TETRAN. 

The  Unit  Tetrahedron 

The  derivation  of  the  exponential  characteristic  method  begins  with  a 
discussion  of  the  unit  tetrahedron  coordinate  system.  This  discussion  will  be 
brief  since  the  face  and  tetrahedra  cell  coordinate  systems,  affine 
transformations,  and  cell  splitting  are  presented  in  detail  in  a  previously 
submitted  paper  (Mathews,  1998).  Consider  the  tetrahedron  in  Figure  3. 


Figure  3.  Tetrahedron  Arbitrarily  Oriented  in  XYZ  Space. 

As  was  pointed  out  by  Brennan  (Brennan,  1996),  it  would  be  an  intractable 

problem  to  develop  a  characteristic  method  using  this  global  cell  coordinate 
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system  because  such  a  method  would  need  to  account  for  any  arbitrary 
streaming  direction  and  cell  orientation.  With  reference  to  Figure  3,  we 
define  an  oblique  Cartesian  (affine)  coordinate  system  with  origin  at  R0  and 
basis  vectors  Et ,  E2 ,  and  E3  for  coordinates  U,  V,  and  W,  respectively.  In  that 

coordinate  system,  the  tetrahedral  cell  is  (by  construction  of  the  coordinate 
system)  always  the  unit  tetrahedron  shown  in  Figure  4. 


Figure  4.  The  Unit  Tetrahedron  in  the  UVW  Coordinate  System. 

Each  cell  in  the  global  mesh  is  mapped  to  its  own  unit  tetrahedron  system. 
All  cell  flux  moments,  source  coefficients,  and  sub-cells  are  defined  based  on 
this  coordinate  system.  After  the  cell  has  been  split  (along  the  streaming 
direction  Q ),  the  inflow  flux  is  passed  from  the  upstream  cell  to  the 
downstream  cell  and  is  properly  transformed  into  the  sub-cell  face  systems. 
This  involves  some  moderate  complexity  and  is  fully  discussed  in  the 
referenced  paper.  In  addition  to  the  cell  coordinate  system,  each  cell  has  four 
faces.  Each  face  has  its  own  coordinate  system  defined  by  two  edge  vectors 
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whose  cross  product  produces  an  outward  normal  vector.  This  face  system  is 
shown  in  Figure  5. 


Figure  5.  The  Face  Coordinate  System. 

With  reference  to  Figure  5,  we  define  an  oblique  Cartesian  (affine)  coordinate 
system  with  origin  at  R£  and  basis  vectors  E[  and  El  for  coordinates  Uf  and 

Vf ,  respectively.  In  this  coordinate  system,  the  tetrahedron  face  is  always 
the  unit  tetrahedron  face  shown  Figure  6. 


Figure  6.  The  Face  Local  UV  Coordinate  System. 

Having  presented  the  needed  coordinate  systems  to  pass  cell  and  face 
flux  moments  within  the  transport  calculation,  we  will  briefly  overview  the 
concept  of  cases  for  cell  splitting. 
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Tetrahedron  Splitting 

The  spatial  quadratures  derived  in  this  work  require  that  each 
tetrahedron  cell  be  split  along  the  streaming  direction  into  a  number  of  sub- 
tetrahedra.  It  is  for  these  sub-tetrahedra  that  the  quadratures  are  derived. 
An  example  cell  split  is  shown  in  Figure  7. 


Figure  7.  Example  of  a  Tetrahedron  Cell  Split. 

Referring  to  Figure  7,  the  splitting  line,  parallel  to  Q ,  enters  the  cell  at 
Rin  =  Ri  and  is  propagated  across  the  cell  to  Rout ,  where  Rout  is  the  point  in 
the  cell  that  will  produce  a  number  of  sub-tetrahedra  that  have  exactly  one 
input  and  one  output  face.  We  term  these  case  1  tetrahedra.  In  the  case 
above,  the  cell  is  split  into  three  sub-tetrahedra  with  corner  nodes 

(R0,R3,Rj,Rout),  ( R3,R2,R1,Rout )>  and  ( R2,R0,R1,Rout )  inXYZ  space  and 
(Xq  ,fj  ,r2 ,r3 ),  (r/,r/,r/,rf),  and  (r/,r/,r/,r/ )  in  each  sub-cell  where  the 

lowercase  refers  to  the  sub-cell  coordinates.  The  node  ordering  is  chosen  to 
define  the  sub-cell  edge  vectors  e1=r1-r0,  e2=r2-r1,  and  e3  =  r3  - r2 , 


20 


which  obey  the  relation  (e1  x  e2)  ■  e3  >  0 .  In  this  case,  the  parent  cell  is  a  case 
-3  tetrahedron  (3  inflow  faces  and  1  outflow  face)  thus  producing  3  sub- 
tetrahedra.  Note  also  that  the  optical  thickness  of  each  sub-cell  is  the  same, 
s  =  a  / ,  where  a  is  the  total  cross  section  for  the  cell  material  and 


1  = 


is  the  length  of  the  splitting  line  (e3  in  each  subcell)  depicted 


in  Figure  7.  Only  six  parent  cell  configurations  for  splitting  an  arbirarily 
oriented  tetrahedron  cell  and  streaming  direction  are  possible.  These  cases 
are  outlined  in  Table  1. 


Table  1.  Possible  Parent  Tetrahedron  Cell  Cases. 


Case 

Sub-cells 

#  of  inflow  faces 
(0 

#  of  outflow  faces 
(/) 

1 

1 

1 

1 

2 

2 

1 

2 

3 

3 

1 

3 

4 

4 

2 

2 

-2 

2 

2 

1 

-3 

3 

3 

1 

Referring  to  Table  1,  the  parent  cell  case  is  determined  by  the  taking  the  dot 
product  of  the  streaming  direction  and  each  outward  face  normal,  n  •  Q .  An 
outflow  occurs  for  n  •  Q  >  0  and  an  inflow  occurs  for  h  •  Q  <  0 .  For  h  •  Q  close 
to  zero  (<  Id12),  no  contribution  is  made  to  the  case  (a  no-flow  face).  Letting  i 
be  the  total  number  of  inflow  faces  and  j  the  total  of  outflow  faces,  the  case  of 
a  tetrahedron  with  respect  to  Q  is  given  as 

Case(iJ)  =  10-11  i  +  i2  -4j  +  5ij.  (3) 
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The  algorithm  used  to  properly  split  a  tetrahedron  into  the  appropriate 
number  of  sub-tetrahedra  is  presented  in  the  aforementioned  paper 
(Mathews,  1998). 

Characteristic  Solution  in  the  Sub-cell  Unit  Tetrahedron 

The  characteristic  solution  of  the  BTE  along  the  streaming  path  Q 
(which  is  parallel  to  the  w  axis —  e3 —  of  the  sub-cell)  is  particularly  elegant 
when  we  operate  on  the  unit  tetrahedron.  Integrating  the  BTE  along  the 
streaming  line,  from  the  point  of  entry  into  the  sub-cell,  at  (u,  v,  0)  (where 
lower  case  u  and  v  are  used  to  indicate  sub-cell  coordinates),  to  the  point  of 
interest,  at  (u,  v,  w),  the  angular  flux  is 

W 

^(u>v,w)  =  \\im(u,v)e~ew+ljS(u,v,w')e~e(~w~w)dw',  (4) 

o 

where  S (u,v,w)  is  the  source,  \\im(u,v)  is  the  inflow  flux,  and  \\i(u,v,w)  is  the 
sub-cell  flux.  Before  proceeding  with  the  derivation  of  the  exponential 
characteristic  method,  we  will  now  briefly  discuss  the  integral  functions 
which  are  central  to  the  EC  method:  the  exponential  moments  functions. 

Moments  Functions 

The  exponential  moment  functions  (Mathews,  1994)  occur  in  deriving 
the  EC  method.  They  are  defined  as 
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%  /  lm-l 

Mn(x1,x2,...,xm)  =  jdti(l-ti)’1  e~Xltl  jdt2e(%1~X2)t2...  jdtm  .  (5) 


They  are  used  to  package  the  numerical  difficulty  of  the  EC  method  and  help 
provide  an  elegant  derivation.  Had  their  pattern  not  been  discerned,  this 
research  would  have  little  hope  of  success.  They  are  related  via  a  multitude  of 
recurrence  relations  and  are  always  positive  (Minor,  1993;  Mathews,  1994; 
Minor,  1995;  Mathews,  1997).  Appendix  A  summarizes  the  exponential 
moment  function  identities  used  in  this  work. 

In  addition  to  the  moment  functions,  various  ratios  of  moment 
functions  are  useful.  One  common  moment  function  ratio  is  the  p  function, 
defined  as: 


At1(xi,x2,---,xIn) 
M0(x1,x2,—  ,xm)' 


(6) 


Note  that  0  <  p(x1,x2,---,xm)  <  1 . 


Lastly,  we  introduce  a  new  moment  function  ratio,  ,  which  is  defined 
as 


^j(Xi,X2,...,Xm) 


M0(x1,x2)...,xm,xj) 

M0(x1,x2,...,xm) 


(7) 


where  the  j  subscript  indicates  that  the  jth  argument  is  repeated  in  the 
numerator  AIq function  argument  list  (but  not  in  the  denominator  A\q).  Note 


23 


also  that  0  <  ^  <  1 .  The  (R]  function  along  with  Identity  6  (Appendix  A)  allow 

for  particularly  elegant  and  robust  formulations  of  the  EC  spatial 
quadrature. 

Key  to  this  research  was  the  development  of  stable  and  accurate 
numerical  routines  for  evaluating  one  through  four  argument  moment  (  M$) 
and  p  functions.  The  one-,  two-,  and  three-argument  function  algorithms 

were  presented  previously  (Mathews,  1997).  The  four-argument  algorithm  is 
presented  in  Appendix  A. 

Derivation  of  the  Exponential  Characteristic  Method 
EC  Source  Distribution 

We  begin  the  derivation  of  the  exponential  characteristic  method  by 
defining  a  volume  moment  operator  for  the  unit  tetrahedron  exactly  as  done 
in  the  linear  characteristic  derivation  (Mathews,  1998): 

l  u  V 

M[g]  =  f  JJJg(E)d3E  =  6jduJdvjdWg(U,v,w).  (8) 

cell  0  0  0 

The  EC  method  uses  an  exponential  source  distribution, 

S(X, Y,Z)  =  exp(a  +  pxX  +  PyY  +  PzZ) ,  that  has  the  known  central  moments 
(initialized  or  known  from  a  previous  iteration),  SA  =  M[S(X,Y,Z)], 

Sx  =  M[(X-X)S(X,Y,Z)],  SY  =  M[(Y- Y)S(X,Y,Z)],  and 

Sz  =  M[(Z  -  Z)S(X.Y,Z)] .  However,  neither  the  moments  nor  their  coefficients 

are  of  any  particular  interest  by  themselves.  Analogous  to  the  LC  derivation, 
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the  moments  produced  by  the  angular  quadrature  can  be  with  respect  to  any 
three  linearly  independent  basis  vectors;  we  choose  the  unit  tetrahedron  edge 
vectors.  Thus,  we  seek  a  function  in  U,  V,  and  W,  that  matches  the  moments 
of  the  source  in  the  oblique  (unit  tetrahedron)  coordinate  system.  The  source 
distribution  then  becomes: 

S(U,  V,  W)  =  exp[As  +  BUU  +  BvV  +  BwW].  (9) 

Operating  on  S(U,  V,W)  with  the  volume  moment  operator,  we  find  the 
following  system  for  the  source  moments: 

SA  =  M[S(U,V,W)]=6exp[AJA)0(Xs,Ys,Zii)  (10) 

Su  =  M[US(U,V,W)]  =  SA  (1  - p(X5,Y8,Zs))  (11) 


Sv  =  M[V S(U,  V ,  W)]  =  Sy  -  SA 


Mq(Xs,Xs.Ys,Z5) 

AVXs,Y„Zs) 


(12) 


and, 


S 


w 


M[WS(U,V,W)]  =  Sv 


Mq(Xs,Ys,Ys,Zs) 

Mo(Xs,Ys,Zs)  ’ 


(13) 


where 


YS=-(BU+BV),  (14) 

Zs  =  -(Bu  +  Bv  +BW). 
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Using  Identity  6  (Appendix  A),  equations  (10)  through  (13)  can  be  written  in 
an  algebraically  efficient  form  as 

SA  =6exp[As]Mo(Xs,Ys,Zs),  (15) 

Sw=SAft3(Xs,Ys,Zs),  (16) 

SV  =  Sw  +  S A  H2  (Xs  , Ys  ,Zs  ),  (17) 

and 

Su=Sv+SAK1(X8,Y.>Zs).  (18) 

Unlike  the  LC  spatial  quadrature,  the  above  system  is  not  directly 
invertible  to  find  the  source  coefficients.  In  order  to  obtain 
As,  Bu ,  Bv,  and,  Bw ,  use  of  a  non-linear  root-solver  is  necessary;  we  use 
Broyden's  method  (Burden,  1993).  Although  the  development  of  an  accurate 
first  guess  algorithm  and  efficient  root-solving  were  key  contributions  of  this 
research,  the  details  of  these  algorithms  have  been  relegated  to  Appendix  B 
in  order  to  preserve  the  flow  of  the  following  derivations. 

Once  the  coefficients  are  found,  direct  affine  transformations  are  used 
to  transform  the  coefficients  from  the  unit  tetrahedron  cell  to  the  various  unit 
tetrahedron  sub-cells  for  use  in  the  spatial  quadrature.  The  reason  for  this  is 
as  follows.  Consider  two  linear  functions,  l(u(x))  and  l'(x),  where 

l(u(x))  =  l'(x)  and  x  is  a  coordinate  vector.  Now,  suppose  that  l(u)  =  a  +  b-u 
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and  u(x)  =  u0  +  J  -(x  -  x0)  where  u0  is  some  initial  value,  J  is  the  Jacobian 
of  the  transformation,  and  x0  is  a  translated  origin  for  the  subcell  coordinate 

system.  Then  we  have  that  l(u(x))  =  a  +  b-u0  +  b-J  • x-b  J  • x0  =a'  +  b'-x 

where  a'  =  a  +  b  -u0  -  b  J  -x0  and  b'  =  b  ■  J  .  In  the  case  of  the  EC  method, 
because  the  source  coefficients  are  related  linearly  in  the  spatial  coordinate 
as  discussed  above,  exponentiating  l(u(x))  and  l'(x)  has  no  effect  on  the 

transformation  of  the  coefficients.  Indeed,  as  long  as  the  coefficients  are 
maintained  in  a  linear  form,  one  can  perform  any  arbitrary  operation,  f ,  on 
both  l(u(x))  and  l'(x)  and  still  perform  a  linear  transformation  of  the 
coefficients  as  presented  above. 

EC  Inflow  Face  Flux  Distribution 

Analogous  to  the  source  equation  development,  we  must  define  an 
operator  to  take  moments  over  an  inflow  face.  In  the  case  of  the  unit  face,  this 
operator  is 

Mf[g]S-VlJg(R)d2R  =  2Jduf  JdVfg(Uf,Vf).  (19) 

A  face  0  0 

As  with  the  source  distribution,  we  assume  an  exponential  distribution 
of  the  inflow  flux  on  the  face, 

4l(Uf  ,Vf )  =  exp(Af  +B*  Uf  +  BjVf ),  (20) 

such  that  the  moments  over  the  face  are 
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=  Mf[^in(Uf  ,Vf )]  =  2exp(Af  )Mo(Xf  ,Yf ), 

(21) 

^=Mf[ufTi(Uf,Vf)]  =  'Fl[l-p(Xf,Yf)], 

(22) 

and  ^  =MffvfT4(Uf,Vf)]-^  yf  M)(xf>Xf,Yf) 

L  J  Mo(Xf,Yf) 

(23) 

where  the  coefficients  are 

Xf  =  -Bfu 

Yf  =-(By  +B^). 

(24) 

Using  Identity  6  (Appendix  A),  equations  (21)  through  (23)  can  be  expressed 

compactly  as 

=  2  exp(Af  )A(o(Xf  ,Yf ), 

(25) 

T4  =>Pi«2(Xf,Yf), 

(26) 

and 


1'6=>F4+4lK1<Xf,Yf).  (27) 

As  with  the  source  system,  the  system  of  equations  for  the  inflow  flux 
moments  is  non-linear  and  requires  a  root-solver  (Broyden's  method)  to 
obtain  the  needed  coefficients. 
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Derivation  of  Sub-Cell  EC  Spatial  Quadrature 

The  above  source  and  inflow  flux  coefficients  are  rotated/translated 
into  the  sub -cell  coordinate  system  using  direct  affine  transformations 
(Mathews,  1998).  This  is  possible  since  the  coefficients  are  linear  and  are  only 
operated  on  by  the  exponential.  Having  obtained  the  coefficients  for  the  sub¬ 
cell,  we  now  define  sub-cell  volume  and  outflow  face  area  moment  operators 
and  the  sub-cell  characteristic  equation  for  the  flux,  \\i  . 

Analogous  to  the  cell  volume  operator,  the  sub-cell  volume  moment 
operator  is  defined  as 

1  U  V 

m[g(r)]  =  6j  duj  doj  dwg(u,v,w).  (28) 

0  0  0 

We  use  lower  case  for  sub-cell  variables  and  upper  case  for  cell  variables.  The 
sub-cell  characteristic  equation  for  the  sub-cell,  obtained  by  integrating  the 
BTE  along  the  streaming  direction  from  ( u ,  v,  0)  to  the  point  ( u ,  v,  w),  is 

U) 

\\i(u,v,w)  =  \\im(u,v)e~ew  +ljs(u,v,w')e~e('w~w')dw',  (29) 

o 

where  the  l,  the  path  length  across  the  sub-cell,  accounts  for  treating  the 
source  distribution  as  a  point  function.  This  is  done  to  keep  the  units  of  the 
source  (particles/cm3/s)  and  flux  (particle-cm/  cm3/s)  consistent  regardless  of 
the  coordinate  system  of  the  independent  variable  for  position  (R,  R,or  r ) 
(Mathews,  1998).  The  source  distribution  is  assumed  to  be  an  exponential, 
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S (u,v,w)  =  exp(as  +  b suu  +  b^u  +  b^u;) .  Similarly,  the  inflow  flux  is  also 
assumed  to  have  an  exponential  distribution  on  the  sub-cell  face, 
yin(u,o)=  exp(a/-  +  b(lu  +  bfvv). 

The  derivation  of  the  sub-cell  flux  moments  for  EC  starts  by 
calculating  the  average  cell  flux, 

^  subcell  =  m[y(vL^W)]  =  ^ubcell  +  vgedlf 

¥subcen  =  6  exp(af)M0(xf,yf  ,z[),  (30) 

't'A.src11  =  6  l  exp(as)  Mo(xf  ,yf  ,zf  ,wf), 

where  the  coefficients  for  the  sub-cell  inflow  face  (obtained  via  affine 
transformation  from  the  global  face)  are 


vf  -  hf 

X1  -  DU> 

yi  =  ~(bu  +  by), 

zi  =-(bu  +bj)  +  s, 


and  the  source  coefficients  are 


xi=-bsu> 

y!  =  -(bsu  +  bsv), 

zi  =-(bsu+bsv)  +  e, 

w?  =-(bsu+bsv  +  b^), 


(31) 


(32) 


where  the  source  and  inflow  flux  coefficients  are  denoted  by  either  an  s  or  f 
superscript,  respectively,  and  s  is  the  cell  optical  thickness. 
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Each  coefficient  above  is  subscripted  with  1  to  differentiate  it  from  the 
coefficients  for  the  outflow  flux.  This  notation  promotes  compactness.  Note 
that  each  sub-cell  will  have  its  own  different  set  of  coefficients. 


As  with  the  average  flux,  the  u-moment  of  the  flux  is  defined  as 


v|Cbce"  =  m[u  y(u,v,w)]  =  VEST"  +<£?“, 
¥^U=vr,L““[l-p(^,y1t,21t)] 

vS?“=vi^“[i-P(x!,y?,z?,w;)] 

=VAuteu[*i(zi,yf,zf,w;)+«2(xf,y?,z;,w;)+ 

^3(xf,yf,zf,wf)  +  ^4(x^yf,zf,wf)], 


(33) 


where  Identity  6  (for  three-  and  four- arguments)  in  Appendix  A  is  used  to 
obtain  the  second  relationship  for  i(/®u^ce11  and  \|/®UgrcceU,  respectively.  Using 
the  formulation  produced  by  this  identity  allows  us  to  see  that 

subcell  subcell  subcell  .  subcell 
VA  >Wu  >Vw  >¥w 


The  derivation  of  the  u-moment  of  the  flux  is  a  little  more  complex 
than  the  average  and  u-moment  derivations.  In  the  previous  cases,  the 
integral  definition  for  the  moment  functions  follows  directly  from  the 
application  of  the  moment  operators.  In  the  case  of  the  u-moment,  integration 
by  parts  is  used  to  evaluate  the  integrals.  Applying  the  moment  operators, 
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where  we  see  that  the  second  equations  for  v|/®^ceUand  i|/®u&^e11  are  obtained 
by  applying  the  previous  definitions  of  ij/®^®11  and  v|/*usbrcceU  and  Identity  6. 

Finally,  the  ie- moment  of  the  cell  flux  is  derived  by  applying  the 
moment  operator  and  integration  by  parts, 


<lCM1  =  m[u>  y(u,v,w)]  = 

„,subcell  _  ,,, subcell  subcell  ,yj  ,y^  ,Zj  ) 

MVin  -Vu.in  ~  V  A, in  T7T1 7 — 77“ 

Mo(xi,yi,Zi) 

subcell  subcell  subceU  ^^(xl>yi>yi >Z1  >wl) 

Vw, src  x  u,src  VA^rc  w  ,  s  s  q 

Mb(x!,yf,zf,w;) 

=  HfAUsrc1I[^3(xl.yi.zl>wl)  +  ^4(xl.y!>zl>Wi)], 


where  again  the  previous  moments  and  the  moment  function  identities  are 
used  to  arrive  at  the  second  form  of  the  equations. 
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Re-arranging  terms  in  order  to  minimize  the  amount  of  arithmetic 
associated  with  the  contributions  to  the  cell  flux  moments  due  to  the  source, 


one  obtains  the  following  equations: 


^subcan  =&i  explaJ.'H/xj.y^zf.V'l), 

vSSf  =  viufe"(R3(x!,y;,Z!,w!)+^(x!,y!,z;,w!)), 

+VA“bA(x;,y!,z!,wf), 

+ vS£n«i(xi,yf,z;>w!). 


(36) 


Similarly,  the  following  equations  are  found  for  the  inflow  flux  contribution  to 


the  sub-cell  flux  moments, 


■ri'm1  =  6  exp(af  jAVxJ.yf.zf), 

<r  =  rtt'u«3<4.yM>. 

subcell  _  subcell  ,  subcell,®  /  f  f  f\ 
Vv.in  —  V w,in  +M/A,in  K2 (X1  >yi  >Z1  h 

subcell  _  subcell  .  subcell cd  „f  „f\ 

Vu.in  =  Vvfin  +  VA,in  ^(Xj  ,Yi  ,Z1). 


(37) 


The  outflow  face  flux  moments  will  now  be  derived.  The  characteristic 


equation  for  the  flux  on  the  outflow  face  is 


\yout(u,v)  =  \\i(u,v,v)  =  \\im(u,v)e  sv  +ljS(u,v,w')  e  w'\  (38) 

o 


Note  that  the  outflow  face  equation  is  v^ace  =v  =  w  (refer  to  Figure  4).  The 


outflow  face  area  moment  operator  is 
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(39) 


1  u 

mowf  [g  (u,  v,  u)]  =  2  J  du  J  du  g(u,  v,  v). 
o  o 


The  inflow  face  area  moment  operator  is  not  explicitly  defined  because  the 
basis  vectors  corresponding  to  the  inflow  face  flux  are  passed  to  it  from  the 
upstream  neighbor  faces  (outflow  faces). 

Starting  with  the  average  flux  on  the  outflow  face, 

^outface  _  ,  ^outface  + 

^outface  _  2  exp(af )  A1o(x^,y|),  (40) 

=  21  fflq)(a.)A%(x5.y|,zi), 
where  the  face  coefficients  are 


4  =  -bj 

y\  =  -(K  +bj)+e 


and  the  source  coefficients  are 


XS2  =  ~K 

yS2=-(b*+b*)  +  e 
Zi  =  -(b*+b‘+b‘). 

Applying  the  moment  operator  with  respect  to  u, 


(41) 


(42) 
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outface 
t  u 


-  m“'[u  y(u,v,v)]  - 


outface 

h 

vSface  =  vX!fce[  1-P(x|,y|)] 

=  vff&C8[^i(xi,yS)  +  ^(x|,y|)], 

v3£“  =  v  AUScce  [1  -  p(x ! ,  y !  > z! )] 

=  H'A.src  [^l(xLyi,z2)  +  ^2(x2.yl>zl)  +  ^3(x2>y!>z!)]> 


(43) 


where  again  Identity  6  (Appendix  A)  is  used  to  arrive  at  the  second  equation 
for  the  u-moment  of  the  outflow  flux. 

Finally,  applying  the  outflow  face  moment  operator  and  integrating  by 
parts,  the  following  equations  for  the  u-moment  of  the  outflow  flux  are 
obtained: 


^outface  =  m1[pf(w)]  = 

f  f  f 

outface  _  outface  ...outface  ^o(x2>x2>y2) 
Vv,in  ~Vu,i  n  “VA,in  - 


src  > 


M»(*j.y|) 


fit”  «2(4.y2), 


outface  _  outface  outface  >M()(X2’X2  ^2  >z2) 
Tu,src  —  ru,src  YA,src  .  .  ,  «  « 

Mo(xs2fys2,zs2) 

=  V AUsreCe  [7?2  (x  2  -  y  2  >  z2  )  +  ^3  (x  1  >  y  1  > z  2  )]  • 


(44) 


Re-arranging  terms  in  order  to  minimize  the  algebraic  operations,  the 


the  outflow  face  flux  moments  due  to  the  source  are 


vS?  =  2/  exp(as)Mo(x|,y|,z|), 

V  v^rc  CC  =  V  AUsrcCe  (^2  (x2  >y!  >z2  )  +  ^3  (x2  >y|  >z2  ))>  (45) 
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Similarly,  the  outflow  face  flux  moments  due  to  the  inflow  flux  are 

VAUinace  =2  exp(af)Mo(x!,y2), 

vS*"  =  vSS,“R2(x|.y|),  (46) 


Sub-Cell  Conservation  Equations 

The  sub-cell  conservation  equations  are  derived  by  applying  the 
volume  moment  operator  to  the  BTE  for  each  moment,  yielding: 


Average  : 

)  +  E'I’a  bc‘M  =  ISA . 

(47) 

w-moment: 

3(Vrface-o+e<ubce,1=;s,,, 

(48) 

u-moment: 

3(v»«.=e-v,in)  +  svjubcell 

(49) 

^-moment: 

3  -  vy1**"  +  yb““  =  ISW. 

(50) 

Note  in  (50)  that  y°utface  =  \j/°utface  while  i^nface  =  0  thus  leading  to  the 

correct  conservation  relation  (the  vj/^bceUterm  is  obtained  as  a  result  of  the 
integration  of  the  BTE).  The  conservation  equations  are  used  in  the  spatial 
quadrature  to  check  that  the  quadrature  is  numerically  accurate.  In  the 
event  that  the  sub-cells  are  not  too  optically  thin,  these  equations  could  be 
used  to  replace  some  of  the  more  expensive  spatial  quadrature  formulas 
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(Brennan,  1996).  However,  this  was  not  done  in  TETRAN  because  it  can  be 
numerically  ill  conditioned. 

Linear  Characteristic  Spatial  Quadrature 

Before  discussing  the  implementation  of  the  EC  method  in  our  code 
system,  we  briefly  summarize  the  linear  characteristic  (LC)  spatial 
quadrature  as  developed  in  (Mathews,  1998).  Its  performance  is  compared  to 
that  of  EC  in  Chapter  4.  The  LC  method  has  been  implemented  in  some  other 
codes  systems  because  it  is  much  less  expensive  computationally  than  the  EC 
method,  so  it  provides  TETRAN  with  a  more  mainstream  transport 
capability.  Adding  LC  requires  modest  extra  effort,  because  all  of  the  cell 
splitting  and  transformation  algorithms  are  the  same  regardless  of  which 
characteristic  spatial  quadrature  is  used. 

LC  was  the  first  quadrature  implemented  into  the  TETRAN  code 
system.  It  assumes  a  linear  source  and  inflow  flux  distribution  as  opposed  to 
the  exponential  assumptions  for  EC: 

S(U,  V,W)  =  As  +  BUU  +  BvV  +  BwW  (51) 

and 

Tin(Uf,Vf)  =  Af  +ByUf  +B^Vf.  (52) 
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Exactly  as  was  done  with  the  EC  source  function,  spatial  moments  are 


taken  of  the  LC  source  distribution  (equation  (51)).  The  resulting 
relationships  for  the  LC  source  moments  are: 

SA  =  M[S(U,V,W)]  =  A,  +i§„  +igv  +igw,  (53) 

4  Z  4 

Su  =  M[US(U,V,W)]  =  f  A.  +  f  Bu  +|gv  +  ±BW,  (54) 

4  5  5  5 

SV  =M[VS(U.V.W)]  =  lAs+|gu+^Bv+^gw,  (55) 

and  Sw=M[WS(U,V,W)]  =  Ia,+|§u+Abv+-Lbw.  (56) 

Unlike  the  EC  method  where  the  source  coefficients  (As , Bu , Bv ,  Bw )  must 

be  obtained  by  root-solving,  the  system  of  equations  above  can  be  directly 
inverted  to  give  the  source  coefficients: 

As  =16Sa  -20Su,  (57) 

Bu  =-20Sa+40Su-  20Sv,  (58) 

Bv  =—  20S|j  +  40  Sy  —  20S^ ,  (59) 

and  Bw  =  —  20  Sy  +40S^y.  (60) 
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Similarly,  the  inflow  face  flux  coefficients  are  obtained  by  taking  the 


appropriate  moments  over  the  tetrahedron  cell's  face  using  (52): 


vi  =Mf['P',(Uf,Vt)]  =  A'  +|b'  +iBj, 


'F£=Mf[uf^(Ut,Vt)]  =  |Af+|Bf  +1b', 


and 


=  Mf  [vfv£(Uf  ,Vf )]  =  ^Af  +^Bfu  +  ±B*. 


(61) 

(62) 

(63) 


As  was  done  for  the  source  moment  equations,  the  inflow  flux  moment 
equations  can  be  directly  inverted  to  solve  for  the  inflow  flux  coefficients: 


Af  =  9Tj[  -12t£  , 

(64) 

B„  =  -12vPj[  +24T£ -1244, 

(65) 

and 

Bj  =-12t£  +  24'¥§. 

(66) 

Clearly,  obtaining  the  required  coefficents  to  solve  the  source  and  inflow  flux 
systems  is  much  less  numerically  intensive  than  that  for  EC.  Thus,  it  should 
be  expected  that  the  LC  method  would  be  computationally  cheaper  than  EC. 


Linear  Characteristic  Sub-Cell  Spatial  Quadrature 
Proceeding  with  the  LC  development,  the  spatial  quadrature  is  found 

analogously  to  EC  by  taking  spatial  moments  of  (29)  (using  the  operator 

defined  in  (28))  with  the  appropriate  source  and  inflow  flux  functional  forms. 
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In  the  case  of  LC,  the  source  is  assumed  to  be  linearly  distributed, 

S (u,v,w)  =  as  +  b*u  +  b*u+  b^w .  The  inflow  flux  function  is  assumed  to  be 


linearly  distributed,  +  b [Lu  +  b{v . 

The  cell  average  flux  moment  is: 


^  subcell  =  =  ^subceU  +  y-jbc* 

^subcell  _  X'000(g)  +  b^  ilL100(e)  +  b^  K0X0(S)), 

y A^src ^  =  6l(as  Ko, 0,0,0 (8)+bu  ^,0,0,0 (s)  +  k*  ^0,1,0, o(£)+ 
k w  [-^0,0, 1,0  (£)  ~  ^o,o,o,i  (£)])> 


(67) 


where  the  special  function,  K,  is  defined  as 

1  tf  ^m-1 

^i1,i2,...,im(£)  =  Jdiijdt2...  J  dtmtil42...^e-8t“.  (68) 

0  0  0 

K  is  only  a  function  only  of  the  cell  optical  thickness,  e .  Because  s  is  the 
same  for  each  subcell,  the  necessary  K  function  values  need  only  be 
calculated  once,  enhancing  efficiency.  With  the  introduction  of  the  K- 
functions,  the  LC  spatial  quadrature  is  expressed  in  an  elegant  form. 
Algorithms  that  accurately  and  efficiently  calculate  the  needed  ^-functions 
for  the  LC  quadrature  are  discussed  in  the  previously  mentioned  paper 
(Mathews,  1998). 

The  u-,  v-,  and  zu-moment  of  the  sub-cell  flux  are  given  as 
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subcell 
r  u 


,  subcell 
,src  9 


=  m[u  v?(u,v,w)]  =  vSe11  +  M'S 
Vw.in  =  i?1>oo(s)  +  bf<  ^2,0,o(£)  +  by  -^1,1,0 (e))’ 

M'u.src  =  6Z(as  -^^0,0,0 (£)  +  bu  -^2,0,0,o(£)  +b*  ^i,i;0,o(£)  + 

bL[^l(0,l,o(£)  “  -^1, 0,0,1 (£)])> 


(69) 


vZnb“U  =  m[in|/(u,i;,i0)]  =  yffi*  +  y", 

VVin  =  6(a^-  K oio (s)  +  b£  ■Z^i>i(o(8)  +  b{  -Ko,2,o(£))> 


subcell 
■  u,src 


6/(as  ^o,i,o,o(s)  +  ^a  Klx0fi(e)  +  b^  ^o,2,o,o(8)  + 


(e)-JKr0,i,0,i(e)]>, 


(70) 


and 


subcell 
t  w 


,  subcell 
>,src  9 


=  m[wy(u,v,w)]  =  v|/^nceU  +  \|/£j 
v"  =  e(a,  ^0,0,1^)  +  ^  KlfiX(z)+h{  tf0,i,i(4 
vL.src11  =  6Z(a s  ■^o,o,i,o(£)+bu  ^i;o(i>0(e)  +  b®  -Ko,i,i,o(£)+ 
bw[-^o,o,2,o(£) -  ^o,o, 1,1  (£)])- 


(71) 


Operating  on  (38)  with  the  outflow  face  area  operator  (39)  using  the 


linear  functions  presented  at  this  beginning  of  this  section,  the  LC  outflow 


face  flux  moments  are 


^outface  =  v  ^ _  ^catfac  +  ^outfte 

rt?"  =2(*fK0,0(e)+b>uK1:0(s)+b(,K0: x  (E>) ,  (72) 

'KaIstc™  =  (as  if000(s)+b^  JC1A0(s)  +  b^,ji?010(e)-  bj)  K0  0  ,(s)); 
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outface 
T  u 


mout[u\y(u,v,v )]  =  y 


outface  ,  outface 
u,  in  ~l"Tii>src  > 


outface 
T  u,in 

outface 
T  u,src 


-  2^a f  if10(s)  +  b„  iiL2,o(e)  +  ^  Kit1"-)) >  (73) 

=  2i  (as  ^i,o,o(s)+bfi  ^2,o,o(8)  +  [K  +  b«j]  ^1,1,0  (s)-  K.  tfl,o,l(8)). 


and 


vo«dace  _  m°“‘[m|/(u,u,a)]  =  y  S?'8  +  <£?”. 

v|C?a“  -  2(a/  J4T0jl(6)+bi  Ku( e)+b'  K02(e)),  (74) 

=  2 1  (as  iST0.l,o(e)+l>J  Kx ls>(z)  +[b“„  +  ba  ]ir0>2>0(s)-  K  JC„,u(e)). 

Although  the  source  and  inflow  flux  moment  coefficients  are  easily 
obtained  for  LC  as  opposed  to  EC,  it  is  not  obvious  from  the  above  moments 
that  the  method  is  numerically  cheaper  than  EC.  However,  what  is  not 
shown  above  is  that  the  various  needed  K- functions  are  related  by 
recurrences  such  that  only  a  handful  of  functions  need  be  calculated.  The 
needed  functions  are  obtained  from  stable  recurrences.  See  (Mathews,  1998) 
for  details.  Additionally,  the  .^-functions  are  functions  of  one  variable 
whereas  the  moments  functions  needed  for  EC  have  two,  three,  or  four 
arguments. 

This  chapter  began  by  introducing  the  notion  of  characteristic 
methods.  Specifically,  the  exponential  characteristic  method  was  derived 
assuming  an  exponential  distribution  of  the  source  and  inflow  flux.  In  order 
to  solve  the  EC  quadrature  for  a  spatial  cell,  non-linear  root-solving  must  be 
done  to  obtain  the  unknown  source  and  flux  parameters.  Following  the 
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discussion  of  the  source  and  inflow  flux,  we  derived  the  EC  spatial 
quadrature  formulas  to  calculate  the  average  and  first  spatial  moments  of  the 
angular  flux  for  a  sub-cell  which  are  obtained  by  splitting  the  parent  cell.  The 
linear  characteristic  method  was  also  briefly  presented  along  with  its  spatial 
quadrature  formulas.  Several  special  functions  were  defined  which  package 
the  numerical  difficulties  of  both  methods  and  allow  for  elegant  derivations. 
With  this  mathematical  foundation,  the  stage  is  set  to  discuss  the 
implementation  of  the  EC  and  LC  quadratures  in  a  radiation  transport  code 
system. 
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Chapter  III:  TETRAN  Implementation 

TETRAN  was  designed  to  be  the  kernel  of  a  parallel,  general  radiation 
transport  code  for  shielding  applications.  In  this  chapter,  we  first  present  the 
programming  model  used  to  achieve  parallel  operation.  The  algorithm  used 
by  TETRAN  to  solve  multigroup,  anisotropic  scattering  problems  is  then 
discussed.  Finally,  the  pseudo-code  for  TETRAN  is  presented  to  familiarize 
the  reader  with  the  code  implementation. 

Parallel  Programming  Model 

The  primary  objective  of  performing  radiation  transport  using  multiple 
processors  is  to  analyze  large  problems  in  a  reasonable  period  of  time.  By 
partitioning  the  problem  space  among  many  processors  operating  on  their 
own  local  data  and  reporting  their  results  as  needed  to  the  other  processors, 
we  hope  to  reduce  the  wall  clock  time  for  getting  accurate  results  for  the 
problem.  We  believe  that  this  is  the  true  measure  of  quadrature  performance. 
Currently,  researchers  compare  method  performance  by  citing  the  time  to 
calculate  the  angular  flux  values  for  a  phase  space  cell  (results  for  one  spatial 
cell  in  one  streaming  direction  for  one  energy  group  and  iteration).  Although 
computational  efficiency  is  an  important  attribute  of  a  spatial  quadrature,  it 
is  not  the  only  criteria  by  which  to  judge  a  method.  Using  phase  space  cell 
performance  criteria  focuses  on  the  computational  efficiency  of  a  quadrature 
without  regard  to  the  methods  accuracy.  By  implementing  parallel 
algorithms,  we  can  ultimately  shift  the  focus  of  quadrature  performance  to  a 
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more  balanced  view  of  efficiency  and  accuracy.  Thus,  the  issues  of 
computational  cost  of  relatively  expensive  spatial  quadratures,  like  EC  and 
LC,  will  become  moot  as  they  are  incorporated  into  parallel  algorithms. 

The  standard  von  Neumann  iteration  on  the  scattering  source  (Lewis 
and  Miller,  1993)  can  be  thought  of  as  an  embarrassingly  parallel  algorithm. 
In  this  algorithm,  we  sweep  through  the  mesh  for  each  streaming  direction, 
Qn ,  calculating  the  needed  flux  moments  (yA ,\|/U,v|/V,andi|/Win  the  case  of 
EC  and  LC)  for  each  cell  in  the  mesh.  After  all  of  the  streaming  directions 
have  been  computed,  numerical  quadrature  is  used  to  calculate  the  scalar 
flux,  <j> ,  for  each  cell  and  check  for  convergence  of  this  value  as  compared  to 
the  last  iteration.  If  the  flux  is  not  converged,  the  scattering  source  is 
updated  using  the  matrix  multiplication,  SA-U’V’W  =  Y  \uA’U’V’W^ ,  .  where 

IjTl  '  IjTI  n  flZfl  5 

nr 

the  subscripts  i  and  n  refer  to  the  cell  and  streaming  direction,  respectively, 

and  T  is  a  material  (i.e.  cell)  dependent  matrix.  T  maps  all  of  the  group  to 
group  and  angle  to  angle  particle  transfers  that  can  occur  and  is  discussed  in 
more  detail  below.  Thus,  the  algorithm  only  requires  communication  between 
processors  to  calculate  the  flux  and  update  the  source  after  an  iteration. 
These  reduction  operations  (dot  products  for  the  numerical  quadrature  and 
matrix  multiplication  for  the  source)  require  expensive  global 
communications  so  that  all  of  the  processors  can  report  their  data  to  each 
other.  Unfortunately,  inter-processor  communication  affects  the  parallel 
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performance  of  an  algorithm  because  the  communication  time  is  generally 
orders  of  magnitude  slower  (Koelbel,  1994)  than  the  compute  speed  of  a 
processor.  Fortunately,  by  pursuing  this  strategy,  we  make  parallel  the 
expensive  computations  required  to  get  the  angular  flux  moments  during  the 
mesh  sweep.  If  the  scalar  flux  and  source  reductions  are  done  efficiently,  the 
algorithm  should  scale  with  the  number  of  processors  attacking  the  problem 
up  to  the  point  where  inter-processor  communications  dominates  the  cost  of 
the  problem. 

In  order  to  implement  the  above  algorithm,  we  chose  to  use  High 
Performance  Fortran  (HPF).  In  HPF,  data  distribution  directives  (!HPF$ 
ALIGN  and  !HPF$  DISTRIBUTE)  are  used  at  a  high  level  of  abstraction  to 
partition  data  among  processors.  Additional  directives  are  used  to  assert  the 
independence  of  loop  iterations  (!HPF$  INDEPENDENT),  thus  partitioning 
computation  among  processors  as  well.  This  coarse  grain  approach  to 
parallelism  shifts  the  work  in  producing  parallel  code  from  the  programmer 
to  the  compiler  and  allows  for  highly  portable  code.  In  contrast  to  the  HPF 
approach  to  parallelism,  a  fine  grain  approach  requires  the  programmer  to 
code  all  of  the  communication  and  data  distribution  requirements  using 
either  the  message  passing  interface  (MPI)  (Snir,  1996)  or  parallel  virtual 
machine  (PVM)  (Geist,  1994)  library  calls.  The  fine  grain  approach  typically 
produces  machine  code  that  is  faster  than  its  HPF  equivalent  given  the 
maturity  level  of  current  HPF  compilers.  However,  fine  grain  codes  also  take 


46 


much  longer  to  write  than  equivalent  HPF  code  and  are  generally  tuned  to 
specific  hardware  architectures.  For  example,  using  HPF  to  parallelize 
TETRAN  took  only  23  total  HPF  directives.  This  was  possible  because 
TETRAN  was  written  with  strict  adherence  to  the  Fortran  90/95  standard. 
HPF  can  be  considered  a  parallel  version  of  Fortran  90.  Thus,  by  using 
Fortran  90  we  were  quickly  able  to  write  a  parallel  radiation  transport  code. 
Additional  discussion  of  the  HPF  language  can  be  found  in  The  High 
Performance  Fortran  Handbook  (Koelbel,  1994). 

Multigroup  Implementation 

The  multigroup  implementation  of  TETRAN  is  a  variation  of  the 
standard  approach  that  is  readily  parallelized.  The  BTE  is  discretized  in 
energy,  denoted  by  the  g  subscript: 

^  +  =^i,n,g’  (75) 

where  i  is  the  cell  index,  n  is  the  streaming  direction  index,  g  is  the  energy 
group  index,  and  a  is  the  total  cross  section.  The  energy  groups  are  arranged 
in  order  of  decreasing  energy  such  that  the  highest  energy  group  is  1.  The 
solution  strategy  is  as  follows.  First,  group  1  is  converged  by  using  the 
standard  von  Neumann  iteration  on  the  scattering  source  (Lewis,  1993).  This 
is  the  within-group  problem  and  the  iterations  are  termed  inner  iterations. 
Having  obtained  the  scalar  flux  for  this  group  and  converged  angular  fluxes, 
the  down  scatter  contribution  of  group  1  to  the  lower  energy  groups  (2 
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through  G)  is  calculated  using  the  scattering  tensor  T  (to  be  discussed 
below).  The  algorithm  then  proceeds  to  group  2  and  repeats  the  process, 
where  group  2  now  has  a  group  1  contribution  to  its  source  particles.  This 
process  is  repeated  until  all  of  the  groups  have  been  converged.  In  the  case 
where  fission  is  present,  there  is  a  possibility  that  lower  energy  groups  will 
be  a  source  of  upscatter  particles.  This  is  repeated  until  the  group  fluxes 
converge  and  is  known  as  the  outer  iteration.  However,  fission  sources  and 
upscatter  are  not  supported  in  TETRAN  making  the  outer  iteration 
unnecessary.  A  block  diagram  describing  this  process  is  shown  in  Figure  8 
and  more  detail  regarding  the  implementation  of  this  algorithm  is  presented 
later  in  this  chapter. 


Figure  8.  Block  diagram  for  iteration  on  the  scattering  source  algorithm. 
Anisotropic  Scattering 

The  traditional  approach  for  modeling  particle  scatter  is  the  spherical 
harmonics  approach,  which  uses  Legendre  moments  of  the  scattering  cross 
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sections  to  model  the  angular  dependence  of  particle  scatter.  The  accurate 
treatment  of  scattering  is  very  important  to  discrete  ordinate  calculations 

because  the  BTE  is  solved  by  decomposing  the  computational  domain  in  Q 
(equation  (2))  and  the  scalar  fluxes  are  calculated  using  numerical 
quadrature.  If  we  use  a  poor  model  or  bad  data,  the  angularly  dependent 
scattering  sources  will  be  wrong  and  the  calculation  will  converge  to  the 
wrong  solution.  The  spherical  harmonic  expansion  of  the  scattering  source  is 

L  l  G 

S,(FA>2:  X  Y/m(6n)X<^(rHg'(r),  (76) 

l-0m=~l  g'-l 

where  S  is  the  scattering  source,  Yim  is  the  conjugate  spherical  harmonic, 

®lgg'  are  the  Legendre  moments  of  the  scattering  cross  section  from  group  g' 

to  group  g ,  and  (j^”.  are  the  coefficients  of  a  spherical  harmonic  expansion  of 
the  group  flux, 

N 

*&<?)  =  Xw(QB0Y/m(QB0y*(r,nn.),  (77) 

7l'  =  l 

where  the  quadrature  approximation  in  (2)  is  used  to  calculate  the  moments. 
Substituting  (77)  into  (76)  and  rearrange  the  summations,  one  arrives  at  the 
following  expression, 

S^(r,Q;i)  =  —>Qn,r),  (78) 

n'  g' 
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where  the  elements  of  the  matrix  T  are 

g'^g 

Tg^g&n'  -> On,r)=w(Qn,) 'Z<ysl>g^g(r)  £  Y^QJY^Q,,,).  (79) 

/-0  m=~l 

Although  T  is  conceptually  a  tensor  that  maps  all  particle  energy  and 
direction  tranfers  in  the  momentum  space,  it  is  used  as  a  matrix  in  the 
following  manner.  Each  material  in  the  problem  has  a  unique  T 'n’->n,g,-*g,mat  > 

where  the  mat  subscript  indicates  the  material  index  where  the  total  number 
of  materials  is  nmat.  During  the  source  update  for  the  within-group  problem, 
r^n’-*n,g^>g,mat  accessed  based  on  the  material  number  of  the  ith  cell  for  the 

current  group,  T n'->n,g'^g,mat(i)  •  This  representation  is  efficient  because  it 

requires  that  only  nmat  different  T  matrices  be  stored  in  memory  during  the 

calculation.  Thus,  it  is  economical  to  replicate  T  across  processors  for  a 
parallel  problem,  further  enhancing  data  locality.  Additionally,  because 
material  regions  tend  to  be  represented  contiguously  in  the  mesh, 
r^n'-+n,g'->g,mat(i)  *s  used  for  several  cells  before  a  new  matrix  is  needed 

because  of  a  change  in  material.  This  reduces  cache  misses  (non-contiguous 
memory  access)  associated  with  reading  in  a  new  T n'->n,g'^>g,mat(i)- 

After  the  within-group  problem  has  converged,  the  resulting  converged 
flux  moments  >  Vu  >  Vv »  an(l  ¥w )  are  used  to  update  the  downscatter 
source  for  the  lower  groups.  This  is  accomplished  as  follows.  The 
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T'n'-*n,g'^>g+i,mat(i) are  calculated  using  (79).  The  source  is  calculated  via  (78) 

and  added  to  contributions  from  other  groups  (accounting  for  the  summation 
in  (76)).  Note  also  that  the  source  terms  are  calculated  and  stored  on  the 
processors  which  access  them  thereby  partitioning  the  problem  in  memory 
among  several  processors.  Thus,  we  are  trading  interprocessor 
communications  for  storage  of  large  problems.  The  above  process  is  done  for 
g  =  g  +  1  to  g  =  G .  Then  the  within-group  problem  is  solved  for  the  next 
group,  and  so  on. 

The  above  approach  is  not  normally  used  since  it  entails  storing  the 
angular  flux  moments  for  every  cell,  which  is  expensive.  For  example,  an  S8 
level-symmetric  quadrature  contains  80  quadrature  directions  which  requires 
(for  LC  and  EC)  storing  80  vj/^ ,  \\ijj ,  v|/y ,  and  v|/^  moments  (320  total)  for 

each  cell  to  calculate  the  scattering  source  moments.  However,  this  algorithm 
is  obviously  parallel  when  the  problem  is  decomposed  along  quadrature 
directions  (Q).  This  methodology  trades  the  benefit  of  decomposing  the 
problem  into  several  sub-problems  running  independently,  each  on  its  own 
processor  (in  its  own  memory  space)  against  the  inter-processor 
communication  required  to  update  the  source  via  the  global  matrix 
multiplication  in  (78).  For  computationally  intensive  spatial  quadratures 
(such  as  LC  and  EC),  we  anticipate  that  this  trade  off  should  result  in  a 
scalable  algorithm  if  the  communications  required  for  the  matrix 

multiplications  are  handled  efficiently. 
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TETRAN  Data  Flow 

TETRAN  requires  geometric,  nuclear,  and  boundary  condition  data  in 
order  to  perform  transport.  This  section  describes  where  the  data  comes  from 
and  how  it  is  assembled  for  TETRAN.  Figure  9  below  shows  the  general  flow 
of  data  for  a  TETRAN  run. 


AVS/Express  format 

Summary  file  current  and  scalar 

flux  data 


Figure  9.  Data  flow  for  TETRAN. 

Pro/Engineer 

Parametric  Technology  Corporation's  Pro/Engineer  ®  18.0  (Pro/E)  was 
used  to  define  the  various  geometries  used  in  this  research.  Pro/E  is  a 
general,  high-end  computer  aided  design  and  manufacturing  (CAD/CAM) 
system,  which  holds  a  large  fraction  of  the  CAD/CAM  market  in  the  United 
States.  Many  national  and  military  laboratories  as  well  as  defense  and 
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energy  contractors  use  it.  Pro/E  excels  in  modeling  large,  complex  assemblies 
quickly  and  has  an  integrated  unstructured  tetrahedron  cell  mesh  generator 
for  use  in  structural  and  thermal  finite  element  analysis. 

Pro/E  was  used  to  perform  several  tasks.  First,  it  was  used  to  develop 
the  individual  parts  that  are  placed  in  a  larger  assembly.  At  the  assembly 
level,  the  Pro/Mesh®  package  was  used  to  add  boundary  conditions  to  each 
assembly  surface.  Additional,  Pro/Mesh®  was  used  to  add  mesh  control  to  the 
assembly  so  that  we  could  perform  mesh  refinement  and  control  how 
individual  parts  were  meshed.  After  meshing,  we  analyzed  statistics 
regarding  cell  aspect  ratios  and  attempt  to  improve  the  mesh.  Pro/Mesh® 
allows  for  mesh  file  output  to  many  finite  element  packages,  such  as 
NASTRAN  (Brennan,  1996).  It  also  has  a  generic  format,  the  FEM  Neutral 
format  (PTC,  1997),  which  is  a  generic  finite  element  mesh  file  giving  all  of 
the  data  needed  by  any  finite  element  code.  We  generated  version  1  FEM 
Neutral  files  to  be  read  by  the  TETRAN  input  module. 

Cross  Section  Data  Format 

The  TETRAN  input  module  uses  TRANSX  2.15,  a  Los  Alamos  National 
Laboratory  cross  section  processing  code,  to  generate  the  problem  specific 
data  libraries  for  TETRAN  (MacFarlane,  1992).  TRANSX  is  a  very  powerful 
cross  section  processing  code  that  can  extract  group-weighted,  self-shielded 
cross  sections  (total  and  Legendre  scattering  moments)  from  a  number  of 
different  libraries.  We  used  it  for  TETRAN  to  mix  cross  sections,  collapse 
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energy  groups,  and  generate  transport  corrected  scattering  moments  in  a 
standard  way.  The  output  file  is  in  the  standard  6E12.5  format,  i.e.,  six 
columns  of  data  per  row,  each  with  a  format  of  E12.5.  Using  TRANSX  allows 
us  to  access  and  process  a  large  number  of  standard  cross  section  libraries  in 
a  well  documented,  straightforward  manner. 

Input  Module 

The  data  for  TETRAN  come  from  a  variety  of  sources.  In  order  to 
reduce  the  complexity  of  TETRAN  and  remove  input  data  processing  from  the 
transport  code,  we  developed  an  input  module  to  read  the  needed  data  and 
put  it  in  a  consistent  format  for  TETRAN.  Thus,  support  for  additional  data 
file  formats  can  be  added  by  revising  or  replacing  the  input  module.  Briefly,  it 
performs  the  following  functions.  For  the  geometry  data,  it  reads  the  FEM 
Neutral  format  file  and  determines  the  number  of  materials  and  boundary 
conditions.  It  then  determines  the  cell  connection  topology,  which  is  an 
important  (but  slow)  step  whereby  the  neighbor  cells  and  faces  for  each  cell  in 
the  mesh  is  determined.  The  connection  topology  is  independent  of  the 
angular  or  spatial  quadrature  chosen  for  the  transport  calculation.  It  is 
needed  in  the  transport  code  to  determine  the  cell  access  order  for  each 
spatial  walk,  which  is  dependent  on  the  angular  quadrature.  Following  the 
generation  of  the  connection  topology,  the  geometry  data  is  written  to  a 
binary  file  to  be  read  by  TETRAN.  Following  the  geometry  data,  the  module 
reads  the  problem  dependent  cross  section  library  created  by  TRANSX  and 
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writes  a  binary  file  for  input  to  TETRAN.  The  input  module  then  reads  any 
user  specified  multi-group  isotropic  source  files  and  writes  this  data  to  a 
binary  file.  Lastly,  for  each  boundary  condition  specified  in  the  mesh  file,  the 
analyst  inputs  the  type  of  boundary  condition  (vacuum  or  incident  current) 
and  the  file  containing  multigroup  incident  current  data.  This  data  is  written 
to  a  binary  file  as  well.  In  addition  to  the  binary  problem  data,  a  header  file  is 
written  which  is  a  text  file  that  tells  TETRAN  several  key  parameters  for  the 
run  such  as  the  desired  spatial  and  angular  quadratures.  These  parameters 
are  changeable  by  the  user  to  allow  for  transport  with  different  angular  and 
spatial  quadratures  for  the  same  mesh  (including  the  connection  topology), 
boundary  conditions,  and  nuclear  data. 

TETRAN  reads  the  above  data  and  performs  the  transport.  Following  a 
successful  run,  the  code  produces  three  output  text  files.  The  first  file  is  a  run 
summary  that  contains  (among  other  things)  the  region  average  scalar  fluxes 
and  boundary  currents  for  the  problem.  The  others  contain  multigroup  scalar 
flux  and  vector  current  data  for  each  cell  in  the  mesh.  This  data,  which  is  in 
AVS/Express  UCD  (Unstructured  Cell  Data)  format  (AVS,  1996),  is  suitable 
for  visualization  with  AVS/Express  (a  scientific  visualization  system). 

Overview  of  TETRAN 

TETRAN  was  written  in  ANSI  standard  Fortran  90,  which  has  been 
the  standard  for  the  Fortran  language  since  1991.  Additionally,  in  some  areas 
of  the  code,  we  have  used  Fortran  95  syntax  (such  as  the  Pure  keyword)  to  aid 
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in  the  parallelization  of  TETRAN.  We  chose  to  adhere  to  the  Fortran90/95 
standard  for  at  least  three  reasons.  First,  Fortran  90/95  contains  many 
intrinsic  functions  such  as  DOT_PRODUCT  and  MATMUL  (general  matrix 
multiplication)  that  allowed  very  readable  and  compact  code  to  be  written, 
aiding  in  the  long-term  maintenance  of  TETRAN.  Second,  Fortran  90/95  has 
adopted  a  superior  model  for  global  data  storage  based  upon  modules, 
allowing  global  data  to  be  passed  throughout  the  code  without  using  the 
traditional  Fortran  77  COMMON  block  (which  was  prone  to  problems). 
Thirdly,  using  Fortran  90/95  allowed  the  quick  addition  of  High  Performance 
Fortran  (HPF)  directives  in  order  to  parallelize  the  code  (Koelbel,  1994). 
Other  reasons  include  Fortran  90/95's  treatment  of  whole  and  masked  array 
operations  and  use  of  explicit  interfaces  for  error  checking.  We  found  the 
book,  FORTRAN  90/95  Explained  (Metcalf,  1996)  and  Digital  Electronic 
Corporation's  (DEC)  Visual  Fortran  5.0  invaluable  in  implementing  TETRAN 
to  the  Fortran  90  /95  standard. 

TETRAN  contains  approximately  9400  lines  of  code  in  87  subroutines. 
The  code  is  intended  to  be  the  kernel  for  a  much  more  capable  transport  code 
system  based  on  the  linear  and  exponential  characteristic  methods  for 
unstructured  tetrahedron  cell  meshes.  This  version  is  capable  of  solving  the 
time -independent  BTE  using  the  standard  source  iteration  technique.  It  is 
capable  of  solving  multigroup  (energy)  problems  assuming  downscatter  only 
(no  fission)  with  general  anisotropic  scattering  based  on  the  spherical 
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harmonics  approach  (discussed  earlier).  Convergence  acceleration,  a 
technique  that  enables  the  source  iteration  algorithm  to  converge  more 
rapidly,  is  not  included  in  this  version.  Additionally,  High  Performance 
Fortran  (HPF)  directives  are  included  for  a  parallel  version  of  the  code.  The 
parallelization  strategy  uses  domain  decomposition  in  the  angular  (Q) 
dimension.  Limited  testing  was  performed  on  the  parallel  version  of  the  code 
due  to  compiler  limitations.  The  code  was  developed  on  a  Hewlett  Packard 
Vectra  XU  2/200  Dual  Pentium  Pro  200  MHz  workstation  using  Digital 
Visual  Fortran,  Version  5.0C.  The  code  was  then  ported  to  the  Aeronautical 
System  Center  (ASC)  Major  Shared  Resource  Center  (MSRC)  IBM  SP  (ASC, 
1998).  We  chose  to  target  the  IBM  SP  because  of  its  large  per  node  memory  (1 
GB),  fast  inter-processor  connection  switch,  and  large  amount  of  per  node 
scratch  disk  (2  GB).  It  compiled  and  ran  serially  under  AIX  4.2  and  the  xlfOO 
compiler,  version  5,  with  minor  modifications.  The  HPF  version  of  the  code 
was  compiled  using  The  Portland  Groups'  PGHPF  2.4  development  compiler. 
All  of  the  data  presented  in  this  report  was  produced  on  one  or  more  nodes  of 
the  ASC  IBM  SP. 

TETRAN  Pseudo-Code 

As  a  concise  and  precise  exposition  of  the  structure  of  the  TETRAN 
program,  we  present  the  following  pseudo-code  outline.  Note  that  HPF 
directives/constructs  are  denoted  by  !HPF$. 
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Begin 


Read  runtime  data  in  problem  data  file:  names  of  binary  input  files, 
spatial  quadrature,  angular  quadrature,  etc. 

Read  binary  data  files. 

Distribute  quadrature  directions  on  to  processor  grid 

! HPF$  distribute  (*, block)  ::  streaming  directions  Qn 
(note:  all  mesh  and  nuclear  data  (cross  sections)  are 
replicated  across  processors) 

a 

Align  transport  arrays  (fluxes  and  sources)  to  Qn  (ensures  other 
computational  arrays  are  on  the  same  processor  to  avoid  unnecessary 
communication) 

! HPF$  align  (*,:)  with  Qn  (*,:)  :: 

V  A  '  V  u  '  Vv  '  Vw  ' 

(each  variable  partitioned  among  processors  along  columns 
of  the  array.) 

! HPF$  Independent:  For  each  quadrature  direction,  fin  ,  determine  the 
cell  access  order  and  case  for  the  spatial  walk. 

For  each  of  the  G  groups  in  the  problem  (g=  1,  2,..,  G;  downscatter 
only)  : 

Read  current  group  cross  sections  and  generate  within  group  T- 
matrix  for  each  material . 

!HPF$  Extrinsic  (HPF_LOCAL) :  Read  group  source  for  each  cell  and 
direction,  (an  extrinsic  subroutine  is  required  to  keep  the 
source  data  i/o  local  to  the  processor) 

Read  current  group  boundary  data  for  each  boundary  cell  and 
direction. 

Do  until  converged  (Inner  Iteration) : 

Set  all  cell  angular  flux  moments  to  zero. 

!HPF$  Independent:  *Meshsweep:  For  each  direction,  Qn  , 
perform  spatial  walk.  (*presented  below)  (the  fluxes  for 
each  quadrature  direction  are  calculated  local  to  their  own 
processor) 

Test  Convergence:  Calculate  scalar  fluxes  and  check 
relative  change  from  last  inner  iteration.  If  not 
converged,  update  within  group  source  for  each  cell  and 
direction  using  T-matrix  for  the  cell's  material.  (1HPF$: 
this  is  where  the  global  communications  occur  using  dot 
products  and  matrix  multiplications.  Note  that  we  let  the 
compiler  determine  the  communication  patterns  for  the 
global  communications . ) 
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<j)(l  rncells)  =  Matmul[ij;(l  :ncells ,  1  :nang)  , weight  (1  :nang)] 
For  i=l  to  ncells 

SA,u,v,w(i.1  :nang)  =  Matraul[v|/AfU|V#w(i,  l:nang),T(l  :nangf  1  :nang,matid  (i))]  + 

®a,u,v,w^  ' 1  :nang) 

Next  cell 

End  do 

Down  Scatter  Source:  Loop  over  current  group  downscatter  cross 
sections,  calculating  T-matrices  and  then  the  current  group 
contribution  to  lower  groups.  !HPF$  Extrinsic  (HPF_LOCAL) :  Update 
the  source  for  each  lower  group  and  store  until  needed,  (again, 
local  routine  is  used  to  keep  sources  local  to  the  processor  they 
are  used  on) 

For  groups  g'=g+l  to  G 

Calculate  T-matrices  for  group  g1  and  materials 
For  each  cell  i 

Sa.u^wUjI  :nang)  =  Matmul[\|/A(UVW(i ,  1  :nang),T(l  :nang,  1  :nang,matid(i))]  + 

SA?u,v,w(i/l:nang) 

Next  i 

I HPF$  Extrinsic (HPF_Local)  Write  to  local  scratch 
disk  for  later  use. 

Next  g' 

Write  Group  Output:  Write  current  group  scalar  fluxes  and 
currents  to  file  for  later  processing. 

Next  Group 

Write  Flux,  Current,  and  Summary  Data 
End 

Of  course,  the  above  is  a  top-level  view  of  the  overall  operation  of  TETRAN. 

The  real  transport  work  is  done  in  the  Meshsweep  loop,  shown  below.  Note 
that  all  of  the  routines  within  Meshsweep  are  Pure  routines,  i.e.  they  produce 
no  side-effects  as  is  required  of  !HPF$  Independent  loops. 


59 


Meshsweep:  !HPF$  Independent:  Do  sequentially  or  in  parallel  for  each 
quadrature  direction,  Qn  . 

Set  face  flux  moments  to  zero  and  initialize  boundary  flux 
values . 


For  each  cell  (in  the  walk  sequence  for  this  direction) : 

Get  cell  index  and  determine  cell  parameters  (case,  source, 
etc. ) 

Determine  cell  split  parameters . 

Get  source  parameters  for  cell: 

LC,  direct  calculation;  EC,  generate  initial 
estimates  then  root  solve 

Determine  cell  upstream  neighbor  faces  and  calculate 
transformation  matrix  to  get  inflow  flux  moments. 

For  each  sub-cell: 


Get  inflow  flux  coefficients  either  directly  for  LC 
or  by  root  solve  for  EC.  Use  direct  affine 
transformation  to  get  coefficients  into  correct  sub¬ 
cell  coordinate  system. 


Get  source  coefficients  by  direct  affine 
transformation  from  parent  coordinate  system. 


Perform  spatial  quadrature  using  correct  sub-cell 
flux  and  source  coefficients  to  get  \|/^  ,  \\fu  ,  i)/^  , 


°ut 

Va 


Wu  ’  and  Vl) 


Vhm 

Check  sub-cell  conservation. 


Next  sub-cell 


Transform/combine  sub-cell  flux  moments  back  to  parent  cell 
UVW  coordinate  system. 

For  each  output  face,  transform/ combine  outflow  flux 
moments  to  parent  global  UV  face  system. 

Next  cell  in  walk 


End  do 


More  detail  regarding  the  stable  evaluations  in  Meshsweep  are  shown  in 
Appendix  C. 

In  this  chapter,  we  discussed  our  strategy  for  the  implementation  of  a 
parallel  radiation  transport  code.  Using  the  spherical  harmonics  source 
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equations  and  the  discrete  ordinates  equation,  we  developed  a  new  approach 
to  calculate  source  updates  that  involves  the  T-matrix.  With  this  approach  to 
the  scattering  source  updates,  the  iteration  on  the  scattering  source 
algorithm  was  implemented  using  Fortran  90/95  in  a  parallel-ready  way 
allowing  the  solution  of  multigroup,  anisotropic  scattering  problems.  High 
Performance  Fortran  (HPF)  directives  were  used  to  decompose  the  needed 
transport  arrays  along  their  angular  dimension,  allowing  the  spatial  walks  in 
all  streaming  directions  to  be  done  in  parallel.  An  overview  of  these 
algorithms  and  the  role  of  the  spatial  quadrature  were  presented  along  with 
the  sources  of  needed  input  data.  With  the  needed  spatial  quadratures  and 
algorithms  developed,  we  now  move  on  to  Chapter  4  where  we  test  TETRAN's 
performance  on  a  variety  of  problems. 
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Chapter  IV:  Testing 

In  this  section,  we  present  the  results  of  our  testing  of  TETRAN  and 
the  exponential  characteristic  method.  Each  test  was  selected  to  demonstrate 
important  features  of  TETRAN's  performance.  The  problems  presented  here 
(except  the  demonstration  of  parallel  code  scaling,  which  was  run  with  2  to  32 
processors)  were  run  on  a  single  node  of  the  ASC  IBM  SP  (known  as  hpc02) 
(ASC,  1998).  Additionally,  the  S8  level- symmetric  angular  quadrature  was 
used  for  all  problems.  The  code  was  compiled  using  optimizations  (-02;  an 
IBM  specific  compiler  flag)  and  tuning  specific  to  the  processor  architecture 
under  the  XLF90  version  5  compiler.  We  compare  the  scalar  flux  and  exiting 
current  results  of  each  problem  with  MCNP  results  using  the  same  input 
data  and  geometry  except  for  the  coarse  mesh,  deep  penetration  problem. 

Note  that  the  current  is  presented  in  units  of  particles/s  through  a  surface. 

We  do  this  because  MCNP  current  data  is  presented  in  these  units.  The  HPF 
version  of  the  code  was  compiled  on  the  ASC  IBM  SP  using  The  Portland 
Groups'  PGHPF  2.4  development  compiler  (an  alpha  version  compiler) 
(Portland  Group,  1996).  The  convergence  tolerance  was  10-6  maximum 
relative  difference  between  the  final  and  next  to  last  iteration  cell  scalar 
fluxes. 

Seven  problems  were  chosen  to  test  specific  aspects  of  TETRAN's 
performance.  Our  first  problem  is  used  for  convergence  testing:  a  uniform 
cube  with  a  uniformly-embedded  isotropic  source.  Our  second  problem  is 
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based  on  problem  1  and  is  used  as  the  basis  for  gathering  of  run-time 
statistics  for  TETRAN's  operation  on  the  IBM  SP.  Problems  3,  4,  and  5 
demonstrate  three  different  aspects  of  robustness  of  TETRAN.  Problem  3  is  a 
uniform  cube  of  material  with  a  uniformly-embedded  isotropic  source 
surrounded  by  a  very  thin  wall  of  shield  material.  This  problem  demonstrates 
TETRAN's  robust  handling  of  very  poorly  shaped  cells.  Problem  4,  a  spherical 
source  within  a  tetrahedron,  demonstrates  TETRAN's  ability  to  perform 
transport  on  meshes  with  curvilinear  and  inclined  plane  surfaces  and 
demonstrates  a  limitation  of  current  mesh  generators.  Problem  5 
demonstrates  the  thick  cell  performance  of  EC  as  compared  to  LC.  This 
problem  is  composed  of  a  uniform  cube  with  a  uniformly  embedded  isotropic 
source  surrounded  on  three  sides  by  an  optically  thick  shield.  Our  sixth 
problem  is  similar  to  one  found  in  the  literature  (Castrianni,  1998),  nested 
cubes  of  water,  iron,  and  water  with  a  source  dissolved  in  the  innermost  cube. 
This  problem  is  used  to  demonstrate  TETRAN's  ability  to  perform 
multigroup,  anisotropic  scattering  problems.  Finally,  problem  7  demonstrates 
scaling  of  our  HPF  parallel  implementation  of  TETRAN  for  the  simple  cube 
problem  used  to  gather  runtime  statistics.  Table  x  summarizes  the  test 
problems  discussed  in  this  chapter  and  the  code  feature  stressed  by  the 
problem. 
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Table  2.  Summary  of  Test  Problems 


Test  Problem 

Feature  Tested 

1 

Convergence  and  convergence  rate 

2 

Run-time  performance 

3 

Poorly  shaped  cell  robustness 

4 

Poorly  shaped  cell  robustness  and  mesh  volume 
conservation 

5 

Robustness  with  optically  thick  cells 

6 

Multigroup,  anisotropic  scattering  perfromance 

7 

Parallel  performance 

Convergence  and  Convergence  Rate 

Our  first  test  problem  demonstrates  the  convergence  and  estimates  the 
convergence  rate  of  the  EC  and  LC  methods  on  simple  meshes.  The  problem 
is  a  1000-cm3  cube  with  a  volume-distributed  isotropically-emitting  source  of 
1.0  particle/(cm3  -  s)  and  vacuum  boundaries.  All  of  the  particles  are  mono- 

energetic.  The  total  cross  section  of  the  cube  material  is  a  =1.0  cm-1  with  an 

isotropic  scattering  cross  section,  cts  =0.5  cm-1 .  TETRAN  was  run  on  seven 
increasingly  finer  meshes  with  nearly  uniformly  shaped  tetrahedra.  The 
characteristics  of  these  meshes  (number  of  cells  and  optical  thickness 
statistics)  are  shown  in  Table  3. 

The  results  produced  by  TETRAN  for  each  mesh  are  listed  in  Table  4 
below  with  the  corresponding  MCNP  and  MCSN  results.  The  MCNP  results 
are  presented  to  show  that  the  method  is  converging  toward  the  expected 
physical  solution  but  are  not  used  to  estimate  convergence  rates  because  they 
contain  no  discretization  errors. 
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Table  3.  Mesh  Parameters  for  Convergence  Rate  Problem. 


Mesh 

emin 

8 

smax 

Emax/ 

/  smin 

1 

6 

8.67230 

17.321 

4.1036 

2 

161 

1.2371 

2.89800 

5.7875 

4.6783 

3 

1292 

0.56301 

1.44680 

3.0510 

5.4191 

4 

4352 

0.37680 

0.96461 

2.0336 

5.397 

5 

10330 

0.26177 

0.72388 

1.5889 

6.0698 

6 

20157 

1.2189 

5.9976 

7 

47794 

0.10087 

0.43431 

1.0050 

9.9633 

In  order  to  determine  the  rate  of  convergence  of  the  methods,  Dr.  Kirk 
Mathews  developed  a  Monte  Carlo  transport  code,  MCSN,  which  transports 
particles  along  the  directions  used  in  discrete  ordinates  angular  quadrature 
sets.  Thus,  the  angular  discretization  error  of  the  discrete  ordinates  angular 
quadratures  is  accounted  for  in  MCSN  allowing  us  to  perform  mesh 
refinement  comparisons  with  the  MCSN  benchmark.  Because  the  MCSN 

benchmark  does  not  contain  spatial  truncation  error  (it  is  continuous  in  R3), 
but  does  contain  the  angular  quadrature  error,  the  relative  error  between  the 
benchmark  and  the  discrete  ordinates  solutions  for  several  meshes  provides 
us  with  our  convergence  rate  graph.  Contrast  this  result  with  the  MCNP 
result,  which  does  not  contain  either  spatial  or  angular  discretization  error 
because  particles  are  sampled  uniformly  in  both  R3  and  Q . 
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Table  4.  Results  for  Problem  1. 


Mesh 

Scalar  Flux 
[particles/cm2  -  s] 

J  ( x=5.0  cm) 
[particles/  si 

1 

1.511537 

40.92211 

2 

1.540445 

38.31644 

3 

1.547338 

37.73447 

4 

1.549319 

37.54896 

5 

1.549909 

37.51102 

6 

1.550281 

37.47919 

7 

1.550566 

37.45456 

MCSN 

1.55091 

(±0.00001) 

37.42417 
(by  conservation) 

MCNP 

1.554800 

(±0.000311) 

36.99300 

(±0.05179) 

Table  4  shows  the  scalar  flux  in  cube  and  the  current  out  of  the  +x  face  at 
x  =  5.0  cm .  Clearly,  both  the  flux  and  current  are  converging  toward  a 
solution  that  is  close  to  the  MCSN  solution. 


To  estimate  the  order  of  convergence  for  the  EC  method,  we  use  the 
analog  to  the  rectangular  cell  case,  assuming  that  average  linear  cell 
dimension,  Ax ,  is  approximated  by  e .  As  the  mesh  is  refined,  e  decreases  by 

a  factor  n  =  ,  where  sx  >  s2.  In  the  thin  cell  limit,  the  mesh  refinement 

/  b2 

should  result  in  a  reduction  in  the  solution  error  by 


n 


p  _ 


(80) 
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where  er  l  and  er  2  are  the  solution  errors  relative  to  a  benchmark  result  for 

the  two  meshes  and  p  is  the  order  of  convergence.  The  absolute  relative  error 
for  the  scalar  flux,  <|) ,  is  defined  as 


,  _  l^il  _  |  4* bench  ~4>i 

^ i  ^  bench 


(81) 


where  <f>  bench  is  the  scalar  flux  predicted  by  MCSN  and  ^  is  the  average 
scalar  flux  in  the  cube  predicted  by  TETRAN  for  the  ith  mesh.  The  relative 
errors  for  both  the  EC  and  LC  quadratures  are  plotted  in  Figure  10  and 
Figure  11  below  versus  the  average  cell  optical  thickness,  s  . 


Figure  10.  Convergence  Graph  for  Exponential  Characteristic  Method. 
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Figure  10  shows  that  the  EC  method  is  a  stable  coarse  mesh  performer  and 
that  it  appears  to  be  at  least  quadratically  convergent  as  compared  to  the  fan 
of  lines  that  represent  ideal  linear,  quadratic,  and  cubic  convergence  with 
respect  to  the  benchmark  and  the  finest  mesh  solution.  However,  the  true 
covergence  rate  is  obscured  because  Pro/Mesh®  cannot  produce  self-similar 
meshes  with  increasing  refinement.  The  EC  and  LC  methods  have  been 
shown  to  be  cubically  convergent  on  such  meshes  for  an  optically  thin 
problem  (Brennan,  1996).  Thus,  we  conclude  that  the  EC  method  as 
implemented  in  TETRAN  is  also  cubically  convergent  and  the  non-unform 
discretization  is  obscuring  this  behavior. 


Figure  11.  Convergence  Graph  for  the  Linear  Characteristic  Method. 
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As  with  EC,  the  convergence  rate  for  LC  is  being  obscured  by  the  lack  of  self¬ 
similarity  in  the  mesh  refinement.  Thus,  it  is  clear  that  future  research  needs 
a  high  quality  mesh  generator  with  the  ability  to  refine  an  unstructured 
mesh  in  a  self-similar  way.  The  LC  results  are  slightly  more  accurate  than 
the  EC  results  compared  to  the  benchmark  (an  effect  seen  in  Brennan's  tests 
as  well).  This  is  because  this  problem  has  a  scattering  source  and  flux  that  is 
concave  down  throughout  the  cube.  The  EC  method  is  attempting  to  model 
this  behavior  with  an  exponential  that  is  concave  up,  which  introduces  more 
approximation  error  than  the  LC  method,  which  uses  a  linear  function.  Thus, 
it  should  be  expected  that  EC  will  perform  less  accurately  than  LC  for  this 
problem.  EC  s  coarse  mesh  capability  is  demonstrated  in  a  later  test  problem. 

Run-time  Profiling 

We  performed  two  levels  of  timing  for  TETRAN.  At  the  top  level  of  the 
code,  we  timed  the  length  of  time  spent  performing  input  and  output^  the 
spatial  walk,  calculating  the  scalar  flux,  testing  convergence,  and  updating 
the  scattering  source.  The  performance  figures  are  for  a  162-cell  mesh  with 
an  S8  angular  quadrature  (80  directions)  which  took  19  iterations  to 

converge.  The  problem  was  run  on  one  node  of  the  ASC  IBM  SP.  The  results 
for  EC  are  listed  in  Table  5.  There  is  no  substantial  difference  between  the 
EC  and  LC  run-time  results  at  this  level  of  profile. 
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Table  5.  TETRAN  Run-Time  Characteristics. 


Action 

%  of  Time  Spent  During  Execution  of  TETRAN 
using  EC  Spatial  Quadrature 

Read  Input  Data 

0.72% 

Generate  Walks  Lists 

0.29% 

Spatial  Walks 

96.51% 

Calculate  Scalar  Flux,  Test 
Convergence,  and  Update 
Source 

0.65% 

Calculate  Currents  and 
Output  Data 

1.39% 

Miscellaneous:  initialize 
arrays,  allocation,  de¬ 
allocation,  other  i/o 

0.45% 

In  addition  to  the  above,  we  compiled  quadrature  specific  data  on  the  cost  to 
perform  transport  using  EC  and  LC  within  TETRAN.  This  data  is  presented 
in  Table  6. 

The  IBM  prof  profiler  (IBM,  1993)  was  used  to  determine  that 
TETRAN  spends  about  45%  of  its  execution  time  evaluating  exponentials  and 
moment  functions  using  EC;  whereas,  the  LC  method  spends  a  negligible 
fraction  of  its  execution  time  evaluating  K  functions.  The  above  performance 
values  (Table  5)  are  for  a  completely  unoptimized  code.  This  is  because  the 
design  philosophy  behind  TETRAN's  development  precluded  the  use  of 
numerical  approximations.  Thus,  TETRAN  serves  as  a  clean  baseline  for 
performance  enhancements  such  as  using  the  conservation  equations  in  the 
spatial  quadrature.  Currently,  the  conservation  equations  are  evaluated  to 
check  the  accuracy  of  the  spatial  quadrature  and  are  not  used  in  the 
transport  calculation  (we  considered  confidence  in  the  results  more  important 
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at  this  stage  than  optimization).  Additionally,  Fortran  90  compilers  are  not 
yet  as  effective  at  optimization  as  Fortran  77  compilers. 

Table  6.  TETRAN  Spatial  Quadrature  Performance. 


Average  Time/Phase  Space  Cell 

(ps/cell-angle-group-iteration) 

(%  of  total) 

Action 

EC 

LC 

Split  cell 

95 

95 

(3%) 

(10%) 

Obtain  source 

537 

34 

coefficients 

(18%) 

(4%) 

Obtain  face  coefficients 

822 

302 

and  translate/rotate  into 

(27%) 

(32%) 

sub -cell 

Spatial  quadrature 

1099 

43 

(36%) 

(5%) 

Re-combine  sub-cell 

200 

200 

moments  into  parent 

(7%) 

(21%) 

cell 

Overhead:  cache  misses, 

260 

260 

page  faults,  resetting 

(9%) 

(28%) 

face  arrays,  timing,  etc. 

Total 

3005 

934 

The  performance  data  in  Table  6  indicate  that  EC  and  LC  are 
moderately  expensive  spatial  quadratures.  However,  note  that  in  the 
previous  effort  (Brennan,  1996),  the  LC  method  was  approximately  10  times 
faster  than  EC  fo'r  the  same  mesh.  We  have  reduced  this  by  a  factor  of  three 
with  anisotropic  scattering  included.  Comparing  the  performance  of  TETRAN 
with  other  codes  at  this  time  is  not  appropriate.  TETRAN  is  not  optimized  at 
all.  The  goal  of  this  research  was  to  develop  EC  and  LC  to  be  robust  and 
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accurate  for  coarse  mesh  problems.  Brennan  pointed  out  (Brennan,  1996) 
that  the  EC  method  pays  for  itself  when  accurate  solutions  are  required  on 
coarse  mesh,  deep  penetration  problems.  Additionally,  our  philosophy  has 
always  been  to  prize  accuracy,  consistency,  and  robustness  over 
computational  speed  since  today's  high  performance  computing  hardware  is 
tomorrow's  antique. 

TETRAN,  like  all  unstructured  mesh  transport  codes,  suffers  from 
poor  cache  performance  with  regard  to  current  high-performance  computer 
(HPC)  architectures.  Currently,  most  HPC  architectures  employ  a  hierarchy 
of  memory  access  based  on  caches.  For  example,  one  node  of  an  IBM  SP 
operates  with  a  135  MHz  Power2  Super  Chip  which  has  a  128  Kbyte  data 
cache,  500  256-byte  cache  lines,  and  256  translation  look-aside  buffers  (TLB). 
The  TLB  is  used  to  translate  between  the  virtual  storage  address  of  a  page  of 
memory  (4096  bytes)  and  it's  current  real  storage  address.  Data  is  accessed 
through  the  cache  lines.  There  is  also  an  instruction  cache,  which  is  32 
Kbytes  in  size.  If  the  needed  data  is  not  located  in  one  of  the  current  cache 
lines,  the  data  from  the  oldest  cache  line  is  flushed,  a  cache  miss  (the  data 
was  not  in  cache)  occurs  and  a  new  cache  line  containing  the  needed  data  is 
read  in,  taking  8  to  12  clock  cycles.  If  the  current  data  is  not  in  the  current 
TLB,  then  the  processor  checks  through  all  of  it's  TLB  to  see  if  it  is  there 
(1Mbyte  of  data).  If  it  is  not,  a  page  fault  occurs  and  the  processor  must  use 
the  operating  system  to  go  to  real  memory  and  access  the  data,  costing  36  to 
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56  clock  cycles.  If  the  data  is  stored  in  arrays  of  unit  stride  and  accessed  in 
this  fashion,  then  cache  misses  will  occur  when  the  data  for  a  cache  line  is 
exhausted  and  another  is  needed.  This  is  the  typical  access  pattern  for  a 
structured  Cartesian-mesh  transport  code.  For  an  unstructured  code,  such  as 
TETRAN,  this  is  a  performance  nightmare.  Data  is  accessed  indirectly  for 
each  cell  and  face.  Thus,  there  is  little  hope  that  even  a  small  to  moderate 
size  problem  will  use  more  than  one  value  from  a  cache  line  or  possibly  a 
page.  Thus,  the  code  expends  a  significant  amount  of  cycle  time  looking 
through  memory  versus  performing  calculations.  We  found  both  the  "P2SC 
Overview"  (Chin,  1996)  and  High  Performance  Computing  (Dowd,  1998)  very 
helpful  in  aiding  our  understanding  of  modern  computer  architectures  and 
the  IBM  Power2  implementation.  Regardless,  this  issue  will  be  significant  for 
unstructured  mesh  transport  for  many  years  to  come. 

Cell  Aspect  Ratio  Robustness 

We  now  examine  TETRAN's  robustness  when  presented  with  a  poorly 
conditioned  mesh  with  many  large  aspect  ratio  cells,  where  the  aspect  ratio  is 
defined  as  the  ratio  of  the  longest  cell  edge  to  the  shortest  distance  of  a  cell 
vertex  to  it's  opposite  face.  These  parameters  are  depicted  in  Figure  12. 
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Figure  12.  Cell  aspect  ratio  parameters. 

To  produce  a  mesh  with  poor  aspect  ratios,  we  used  a  solid  cube  of 
material  (1  cm3)  surrounded  on  all  sides  by  a  thin  layer  of  0.01  cm  thickness. 
Two  examples  of  the  resultant  meshes  are  shown  in  Figure  13. 


Figure  13.  Meshes  for  Robustness  Test. 


This  problem  is  chosen  since  it  serves  as  the  basis  of  a  more  complex 
assembly  of  such  thin  walled  objects.  It  is  also  easy  to  model  in  MCNP  so  that 
we  can  use  the  Monte  Carlo  solution  as  a  reference  value.  Lastly,  the  model  is 
representative  in  proportion  to  the  thin  walls  that  typically  surround  a 
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satellite  (100  mils  about  a  1  ft3  box  ->0.01  cm  shield  about  1  cm3  box) 
(Hilland,  1998),  thus  implying  military  as  well  as  civilian  applicability. 

The  test  problem  for  this  demonstration  is  a  1  cm  x  1  cm  x  1  cm  cube 
with  a  constant  isotropic  source  (monoenergetic)  of  1.0  particle/s.  It  has  a 

total  cross  section  of  a =1.0  cm-1  and  a  scattering  cross  section  of 
crs  =0.5  cm-1 .  The  shield  region  is  0.01  cm  thick  and  surrounds  the  cube.  The 
shield  has  a  total  cross  section  of  c  =  2.0  cm-1  and  a  scattering  cross  section  of 
os  =1.0  cm-1 . 

This  test  case  was  designed  to  demonstrate  the  robustness  of  both  the 
EC  and  LC  methods.  The  maximum  cell  aspect  ratio  for  the  coarse  mesh 
problem  is  144:  the  aspect  ratio  for  80%  of  the  cells  in  this  mesh.  The  fine 
mesh  has  cell  aspect  ratios  distributed  as  shown  in  Figure  14.  Clearly, 
refining  the  mesh  reduces  the  occurrence  of  large  aspect  ratio  cells.  Note  that 
the  aspect  ratio  for  an  equilateral  tetrahedron  is  1.2247  indicating  that  none 
of  the  meshes  for  this  problem  contained  ideally  shaped  cells.  The  meshes  for 
this  problem  (Figure  13)  also  exhibit  a  problem  with  current  mesh 
generators.  Notice  in  the  fine  mesh  on  the  right  that  the  corners  are  meshed 
preferentially  versus  the  interior  of  the  faces.  This  is  presumably  because  the 
finite  element  mesh  generator  is  attempting  to  produce  a  good  mesh  for  a 
thermal  or  structural  response  code.  We  do  not  need  this  type  of  mesh  for 
particle  transport.  Instead,  we  need  a  uniform  mesh  of  well-shaped 
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tetrahedra.  This  is  an  instance  where  a  mesh  generator  optimized  for 
radiation  transport  would  be  useful. 


Figure  14.  Distribution  of  Cell  Aspect  Ratios  for  Fine  Mesh. 

The  coarse  mesh  has  the  minimum  number  of  cells  that  Pro/Mesh 
would  use  for  this  problem.  A  finer  mesh  than  the  4902-cell  mesh  presented 
could  not  be  generated,  because  the  number  of  cells  required  grows 
enormously  with  increasing  mesh  refinement.  The  cell  optical  thicknesses  for 

the  two  meshes  analyzed  are  presented  in  Table  7.  Additionally,  £ma y  is 

/  smin 

shown  because  it  reflects  the  capability  of  the  mesh  generator  to  produce 
uniform  meshes. 
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Table  7.  Cell  Optical  Thicknesses  for  Cell  Aspect  Ratio  Problem. 


Region 

Nee  11s 

smin 

s 

emax 

^max : / 

/  smin 

Source 

6 

0.42209 

0.86392 

1.7321 

4.10 

1671 

0.00108 

0.07965 

0.51769 

479 

Shield 

36 

0.01026 

0.02651 

0.17321 

16.9 

3231 

0.00264 

0.03736 

0.09165 

34.7 

Cells  with  large  aspect  ratios  present  difficult  problems  numerically 
for  EC  because  the  source  and  inflow  flux  coefficients  for  these  cells  are 
generally  large,  negative,  and  close  together.  This  behavior  can  cause 
overflows  if  not  handled  properly  (see  Appendix  C).  Badly  shaped  cells  also 
cause  problems  with  the  cell  Jacobian  evaluation  and  rotation  matrices. 
Thus,  this  test  case  to  a  lesser  extent  shows  the  robustness  of  the  LC  method 
as  well.  The  successful  evaluation  of  this  problem  verifies  that  we  have 
achieved  a  level  of  robustness  that  was  not  available  in  the  earlier  work 
(Brennan,  1996).  Note  that  this  problem  does  not  test  the  coarse  mesh  ability 
of  either  EC  or  LC  as  evidenced  by  the  relatively  small  optical  thicknesses  in 
Table  7. 

We  demonstrate  the  performance  of  the  method  by  comparing  the  LC 
and  EC  results  with  those  obtained  from  MCNP.  This  data  is  shown  in  Table 
8.  This  data  shows  that  the  LC  and  EC  methods  are  converging  toward  a 
discrete  ordinates  solution  that  is  generally  in  agreement  with  the  MCNP 
prediction.  Thus,  TETRAN  provides  accurate  results  even  for  poor  meshes. 
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Table  8.  Results  for  Large  Cell  Aspect  Ratio  Problem. 


Method 

Ncells 

Ncells 

^Shield 

^Cube 

J+ 

(Shield) 

(Cube) 

H/(cm2  -s)l 

n/(cm2  -s)l 

<m 

MCNP 

0.231440 

0.403520 

0.130500 

(±0.000069) 

(±0.000081) 

(±0.000104) 

42  Cells 

36 

6 

LC 

0.239455 

0.385757 

0.131714 

EC 

0.240027 

0.384161 

0.131925 

192  Cells 

144 

48 

LC 

0.234940 

0.398565 

0.131053 

EC 

0.235534 

0.397051 

0.131294 

4902  Cells 

3231 

1671 

LC 

0.234192 

0.400116 

0.130974 

EC 

0.234290 

0.399896 

0.130997 

This  problem  demonstrated  the  difficulty  that  unstructured 
tetrahedral  mesh  generators  can  have  with  thin  regions  in  a  mesh.  Clearly, 
building  robust  algorithms  to  handle  poor  meshes  is  a  way  to  overcome  this 
problem.  Another  approach,  which  was  initially  pursued  in  this  research,  is 
to  develop  a  spatial  quadrature  that  treats  thin  regions  in  a  problem  with  a 
new  spatial  quadrature.  This  approach  is  discussed  in  Appendix  D  and  is 
called  the  surface  cell  algorithm.  In  this  algorithm,  the  thin  regions  are 
collapsed  to  a  surface  that  are  treated  as  computational  cells  with  their  own 
spatial  quadrature.  Thus,  a  mesh  would  have  both  surface  cells  and 
tetrahedra  cells:  a  mixed  mesh.  This  research  was  abandoned  because  there 
were  no  mesh  generators  available  that  produced  the  needed  mesh.  Should 
such  a  mesh  generator  become  available  it  would  enable  the  implementation 
of  this  approach. 
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Curvilinear  and  Inclined-Plane  Surfaces 

The  problem  presented  in  this  section  demonstrates  TETRAN's  ability 
to  perform  transport  on  complex,  multi-region  geometries  with  curvilinear 
and  inclined-plane  surfaces.  This  problem  demonstrates  the  robustness  of  the 
algorithms  and  an  issue  regarding  mesh  generation  of  sources  in  a  complex 
configuration.  We  initially  believed  Brennan's  development  able  to  handle 
such  problems.  However,  it  was  found  that  the  code  was  generally  not  robust 
enough  to  handle  the  few  poorly  shaped  tetrahedra  that  are  produced  in  such 
a  configuration.  Consequently,  Brennan's  code  typically  failed  with  overflow 
errors  (see  Appendix  C)  during  execution  for  simple  multi-region  problems. 
EC  was  made  robust  by  using  the  implementation  presented  in  Chapter  2. 

This  problem,  a  spherical  source  within  a  tetrahedron,  was  chosen  to 
demonstrate  the  three-dimensional  capability  of  TETRAN.  Indeed,  only  one 
other  discrete  ordinates  code,  ATILLA,  is  even  capable  of  performing 
transport  on  this  configuration  with  such  a  high  degree  of  fidelity.  MCNP  is 
also  capable  of  performing  transport  on  the  configuration  and  was  used  as 
the  benchmark  for  this  problem.  The  geometric  configuration  is  as  follows. 

For  the  tetrahedron,  the  corners  are  located  in  (X,Y,Z)  space  at  the  following 
nodes:  node  0  (5.0,  3.0,  -3.0),  node  1  (0.0,  0.0,  5.15),  node  2  (-5.0,  3.0,  -3.0),  and 
node  3  (0.0,  -5.66025,  -3.0).  The  sphere  is  centered  at  the  origin  and  has  a 
radius  of  1.336505  cm  (a  volume  of  10.0  cm3).  The  total  cross  section  for  both 
regions  is  0.75  cm1  and  the  scattering  cross  section  is  0.5  cm1.  These  cross 
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sections  were  chosen  because  they  were  used  by  Lathrop  (Lathrop,  1971)  in 
his  square-in-square  problem.  We  envision  this  problem  as  a  three- 
dimensional  version  of  the  square -in-square  problem  (actually,  a  cube-in-cube 
is  the  three-dimensional  variant  but  we  wanted  a  more  challenging 
geometry).  A  monoenergetic,  isotropic  source  of  strength  1.0  particles/(cm3-s) 
is  located  within  the  sphere. 

We  examined  this  transport  problem  using  three  meshes:  a  coarse 
mesh  with  138  cells,  a  finer  mesh  with  1209  cells,  and  the  finest  mesh  with 
10075  cells.  Table  9  lists  the  region  volumes  and  optical  thicknesses  for  each 
region  and  mesh. 

It  can  be  seen  from  the  below  data  that  the  source  volume  for  each 
mesh  is  not  being  conserved  nor  is  that  of  the  tet.  The  actual  source  volume  is 
10.0  cm3.  Although  the  overall  volume  of  the  problem  is  conserved,  if  the 
source  volume  is  incorrect,  the  discrete  ordinates  solution  will  be  incorrect 
because  the  source  will  be  in  error.  In  addition  to  not  having  the  correct 
volume,  none  of  the  meshes  have  the  correct  source  shape  (a  sphere)  although 
they  are  approaching  it  with  each  mesh  refinement.  This  is  shown  in  the 
figures  below,  where  we  present  the  coarse  and  finest  meshes. 
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Table  9.  Mesh  Optical  Thicknesses  for  Curvilinear  and  Inclined-Planes 

Problem. 


Region 

Volume 

(cm3) 

Ncells 

^min 

8 

^  max 

^max : / 

/ 6  min 

Source 

2.4028 

8 

0.44200 

0.76068 

1.5635 

3.54 

8.4789 

186 

0.41175 

9.44 

9.4916 

938 

0.24471 

itlifrVJMif 

9.70 

10.0 

Actual 

N/A 

N/A 

N/A 

N/A 

HUB 

115.23 

130 

0.17931 

1.0240 

2.7929 

15.58 

109.15 

1023 

0.08732 

0.50121 

2.0958 

24.00 

108.14 

9137 

0.05553 

0.26430 

0.62608 

11.27 

107.6351 

Actual 

N/A 

N/A 

N/A 

N/A 

Figure  15.  Coarsest  Mesh  of  Sphere  in  Tetrahedron  (138  cells). 

Clearly,  the  mesh  in  Figure  15  is  completely  unrealistic  for  the  embedded 
spherical  source.  Contrast  this  with  the  finest  mesh,  which  is  significantly 
different  than  the  sphere  it  is  modelling  with  regard  to  the  volumes 
(9.4916  cm3  vs.  10.0  cm3)  (Figure  16). 
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Figure  16.  Finest  Mesh  Sphere  in  Tetrahedron  Mesh  (10075  cells). 

We  cannot  rely  on  the  mesh  generator  to  conserve  volumes  automatically. 
This  feature  is  critical  to  nuclear  applications  because  one  often  deals  with 
volumetric  sources.  In  the  case  of  this  problem,  even  finer  meshes  are 
required  to  get  the  source  right,  although  we  could  multiply  the  region  fluxes 
by  a  correction  factor  to  account  for  the  error  in  the  source  volume  in  the 
mesh.  Such  an  approach  is  plausible  for  a  simple  problem  but  nearly 
impossible  for  a  complicated  source  configuration.  A  volume-conserving  mesh 
generator  is  needed.  The  results  of  the  transport  calculation  (using  EC)  for 
each  mesh  are  presented  below.  Note  that  each  face  is  numbered  according  to 
the  node  that  does  not  lie  on  it,  e.g.  face  0  is  comprised  of  nodes  1,  2,  and  3, 
and  so  on. 
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Table  10.  Results  of  Transport  Calculation  Using  EC  for  Test  Problem  2. 


/  particles  x 

*  2  ' 
cm  -s 

Current  Out  of  Face  :  (particles/s) 

Ncells 

^Source 

^Tet 

0 

1 

2 

3 

138 

0.625647 

(2.60382) 

0.033598 

(0.139830) 

0.364600 

(1.517396) 

0.100795 

(0.419490) 

0.363604 

(1.513250) 

0.230081 

(0.957552) 

1209 

0.979378 

(1.155076) 

0.096365 

(0.113653) 

1.168384 

(1.377990) 

0.377978 

(0.445786) 

1.152992 

(1.359837) 

1.073805 

(1.266444) 

10075 

1.017173 

(1.071656) 

0.104814 

(0.110429) 

1.308251 

(1.378325) 

0.424641 

(0.447386) 

1.306480 

(1.376459) 

1.204846 

(1.269381) 

MCNP 

1.034600 

(±0.000207) 

0.108750 

(±0.000022) 

1.382700 

(±0.000968) 

0.446760 

(±0.000536) 

1.383300 

(±0.000968) 

1.274700 

(±0.000892) 

Notice  that  if  we  multiply  the  scalar  fluxes  and  currents  by  the  source  volume 
ratio  (10.0/source  volume) ,  the  transport  solutions  become  remarkably  more 
accurate  with  respect  to  the  Monte  Carlo  solution.  The  corrected  values, 
shown  in  parenthesis,  are  tabulated  in  Table  10.  Note  that  correcting  the 
source  flux  produces  a  worse  estimate  of  the  scalar  flux  for  the  coarse  mesh 
than  for  the  other  meshes.  This  is  because  the  source  region  in  the  coarse 
mesh  doesn't  have  the  shape  of  a  sphere.  Thus,  it  should  be  expected  that  the 
volume  correction  will  not  improve  the  flux  prediction  in  this  case. 


Clearly,  correcting  for  the  missing  source  due  to  the  mesh  improves  the 
accuracy  of  the  discrete  ordinates  solution  with  respect  to  the  MCNP 
solution.  Note,  however,  that  in  most  complex  source  configurations  it  will 
not  be  easy  to  perform  such  a  correction.  The  better  approach  for  this  type  of 
problem  is  to  mesh  the  problem  correctly  such  that  volumes  are  conserved. 
This  eliminates  any  need  for  such  corrections  and  frees  us  from  worrying 


about  such  problems. 
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The  above  table  shows  that  the  source  correction  works  remarkably 
well  for  some  aspects  of  the  problem.  Rather  than  under-predicting  the  flux 
and  currents,  we  over-predict  and  converge  toward  the  MCNP  solution.  This 
behavior  is  far  more  appealing  for  shielding  calculations  than  the  under¬ 
predicting  behavior  that  occurs  when  the  wrong  source  is  used.  This  problem 
shows  that  work  needs  to  be  done  in  improving  mesh  generation  for 
applications  that  need  good  control  over  the  mesh  volumes.  Without  such  an 
ability,  future  use  of  unstructured  mesh  codes  will  be  limited  since  the 
primary  benefit  of  such  methods  is  to  model  curved  geometries.  The  issue  of 
volume  conservation,  coupled  with  the  finite  element  artifacts  discussed  with 
respect  to  the  last  test,  more  than  drive  home  the  need  for  a  transport  specific 
mesh  generator. 

With  evidence  that  TETRAN  is  robust  and  accurate  (given  the 
characteristics  of  the  meshes  used)  for  complicated  multi-region  problems,  we 
proceed  to  our  next  test,  a  problem  that  demonstrates  TETRAN  ability  to 
solve  coarse  mesh,  deep  penetration  problems. 

Coarse  Mesh  with  Optically  Thick  Cells 

Our  last  test  of  TETRAN's  robustness  centers  on  its  ability  to  solve 
problems  with  cells  that  are  optically  thick.  EC  has  been  shown  superior  to 
LC  for  these  types  of  problems.  The  LC  method,  a  coarse  mesh  linear  method, 
is  reliably  accurate  for  cell  optical  thicknesses  on  the  order  of  2.0  or  less.  The 
difference  methods  (diamond  difference  et  al.)  require  much  finer  meshes.  We 
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demonstrate  here  EC's  performance  for  cells  with  optical  thicknesses  up  to 
34.  Because  this  problem  is  a  stressing  deep  penetration  problem,  we  are 
unable  to  use  MCNP  to  obtain  a  bench  mark  solution  to  the  problem.  The 
thickness  of  the  shield  is  on  the  order  of  20  mean  free  paths.  Rather  than 
experiment  with  variance  reduction  techniques,  we  use  the  most  converged 
EC  solution  as  our  benchmark  for  this  problem. 

This  problem  is  composed  of  two  regions.  The  inner  region  is  a 
homogeneous  cube  (dimensions  10  cm  x  10  cm  x  10  cm)  which  contains  a 
homogeneously  distributed  source  of  monoenergetic  particles  emitted 
isotropically  with  a  strength  of  1000  particles/cm^.  The  total  and  scattering 
cross  sections  of  the  source  cube  material  are  c  =  0.25  cm-1  and 

CTs  =  0.0833cm-1 ,  respectively.  The  outer  region  is  a  20  cm  x  20  cm  x  20  cm 

cube  with  a  =2.0  cm  Jand  a  =0.5 cm-1 .  The  source  material  is  located  in  one 

corner  of  the  larger  cube,  as  is  shown  in  Figure  17.  We  chose  cubes  to  ensure 
that  the  volume  of  the  source  region  was  represented  exactly. 
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Figure  17.  Geometry  for  Coarse  Mesh  Problem. 

The  region  scalar  fluxes  and  exiting  face  currents  were  calculated  for 
four  meshes  using  EC  and  LC.  The  coarsest  mesh  had  48  cells  while  the  finer 
meshes  had  348,  3063,  and  20165  cells,  respectively.  The  mesh  optical 
characteristics  are  presented  in  Table  11. 


Table  11.  Mesh  Optical  Thicknesses  for  Coarse  Mesh  Problem 


Region 

Ncells  in  region 

emin 

s 

emax 

£max/ 

/  smin 

Source: 

6 

1.0552 

2.1763 

4.3301 

4.10 

48 

0.48656 

1.0825 

2.1663 

4.45 

381 

0.22113 

0.54329 

1.0836 

4.90 

2047 

0.11241 

0.30941 

0.69405 

6.17 

Shield: 

42 

8.0038 

17.376 

34.641 

4.33 

336 

3.5912 

8.6710 

17.321 

4.82 

2682 

1.7139 

4.3362 

9.0774 

5.30 

18118 

0.56092 

2.2600 

5.1109 

9.11 

Table  11  shows  that  the  meshes  examined  contain  extremely  thick  cells  in 
the  case  of  the  coarse  mesh  (up  to  e~34)  to  only  moderately  thick  cells  in  the 

case  of  the  fine  mesh  (but  still  thick  enough  to  reduce  LC's  accuracy).  The 
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source  region  is  much  less  optically  thick  than  the  shield  region.  The  results 
for  the  region  average  scalar  fluxes  are  listed  in  Table  12. 


Table  12.  Region  Scalar  Fluxes  for  Coarse  Mesh,  Deep  Penetration  Problem 


^  source 

^  shield 

,  particles , 

^  2  ' 

,  particles . 

V  9  J 

cm 

-s 

cm 

-s 

Mesh: 

EC 

LC 

EC 

LC 

coarse 

2.8239E+03 

2.8000E+03 

2.3757E+01 

2.4895E+01 

(48  cells) 

less  coarse 

2.8993E+03 

2.8816E+03 

2.3167E+01 

2.3709E+01 

(384  cells) 

finer 

2.9185E+03 

2.9101E+03 

2.3041E+01 

2.3322E+01 

(3063  cells) 

finest 

2.9229E+03 

2.9199E+03 

2.3005E+01 

2.3102E+01 

(20165  cells) 


From  the  above  table,  it  can  be  seen  that  both  EC  and  LC  are  converging  to 
the  same  result  for  the  finest  mesh.  The  EC  method  is  more  accurate  than  LC 


for  all  of  the  meshes.  The  relative  errors  with  respect  to  the  finest  mesh  EC 
scalar  flux  for  both  regions  are  plotted  in  Figure  18  and  Figure  19. 


Figure  18  shows  that  the  EC  method  is  more  accurate  than  the  LC 
method  even  in  the  optically  thin  source  region.  However,  LC  is  still  doing  a 
very  good  job  of  getting  the  scalar  flux  in  this  region  because  the  cells  are  not 
optically  thick. 
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Figure  18.  Relative  Errors  in  Source  Region  Scalar  Flux  with  Respect  to  EC 
Fine  Mesh  Benchmark  for  Deep  Penetration  Problem. 


8 


Figure  19.  Relative  Errors  in  Shield  Region  Scalar  Flux  with  Respect  to  EC 
Fine  Mesh  Benchmark  for  Deep  Penetration  Problem. 
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In  the  case  of  the  shield  region  (Figure  19),  the  LC  method  is  significantly  in 
error  with  respect  to  EC.  This  should  be  expected,  however,  because  the 
region  is  optically  very  thick.  LC  has  little  hope  of  matching  the  highly 
exponential  behavior  of  the  flux  in  outer  reaches  of  the  shield  whereas  EC  is 
doing  an  excellent  job  in  this  regard  as  evidenced  by  the  scalar  flux  and 
currents. 

Lastly,  the  exiting  face  currents  for  this  problem  are  listed  in  Table  13. 

Table  13.  Exiting  Face  Currents  for  Deep  Penetration  Problem. 


J+x  J+y  J+z 

(particles/s)  (particles/s)  (particles/s) 


Mesh 

EC 

LC 

EC 

LC 

EC 

LC 

48 

1.1105E-04 

-1.8645E+02 

9.1345E-05 

-1.8316E+02 

9.1895E-05 

-1.8531E+02 

384 

4.1948E-05 

7.2448E+00 

4.2087E-05 

3.9853E+00 

4.064 IE-05 

4.8650E+00 

3063 

3.5524E-05 

-6.6754E-04 

3.5078E-05 

1.2391E-03 

3.5282E-05 

4.0023E-04 

20165 

3.3825E-05 

1.3047E-05 

3.3890E-05 

1.2830E-05 

3.3825E-05 

1.3599E-05 

J- 

X 

J-y 

J-z 

(p  articles/s) 

(particles/s) 

(particle  s/s) 

Mesh 

EC 

LC 

EC 

LC 

EC 

LC 

48 

9.3634E+04 

9.1195E+04 

9.3127E+04 

9.0643E+04 

9.3127E+04 

9.0644E+04 

384 

9.1186E+04 

9.0434E+04 

9.1248E+04 

9.0092E+04 

9.1091E+04 

9.0242E+04 

3063 

9.0545E+04 

8.9966E+04 

9.0554E+04 

9.0071E+04 

9.0542E+04 

9.0057E+04 

20165 

9.0434E+04 

9.0280E+04 

9.0429E+04 

9.0242E+04 

9.0427E+04 

9.0256E+04 

From  Table  13,  we  see  that  the  exiting  face  current  prediction  for  the  far 
faces  of  the  shield  (+x,  +y,  and  +z)  show  the  deep  penetration  performance  of 
EC  versus  LC.  For  the  coarsest  mesh,  LC  produces  returns  negative  currents 
through  all  of  the  shield  faces.  This  is  nonsense  and  we  expect  it  because  of 
the  coarseness  of  the  mesh.  The  EC  calculation  predicts  currents  that  are 
only  a  factor  of  three  higher  than  the  most  converged  solution.  Note  that  the 
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symmetry  of  the  problem  implies  that  J+x  =  J+y  =  J+z  (and  J_x  =  J_y  =  J_z). 

This  behavior  is  clearly  occuring  as  the  mesh  is  refined.  The  LC  current 
prediction  is  obviously  completely  in  error  until  the  finest  mesh  and  then  it  is 
still  considerably  off  because  the  mesh  in  the  shield  is  still  optically  thick 
(s  =  2.26)  for  LC.  Further  mesh  refinement  should  bring  the  LC  result  in  line 
with  EC's.  The  relative  errors  (with  respect  to  the  most  converged  EC 
solution)  in  J+x  are  plotted  in  Figure  20  versus  e  for  the  shield  as  it  is  the 
shield's  optical  characteristics  that  are  effecting  the  computation  of  J+x  . 


•  EC 
O  LC 


Figure  20.  Relative  Error  with  Respect  to  Fine  Mesh  EC  Benchmark  for  J+x 

using  EC  and  LC. 

Notice  in  Figure  20  that  the  current  predicted  by  the  384  cell  mesh  using  EC 
is  more  accurate  than  the  LC  results  on  the  finest  mesh  (20165),  a  factor  of 
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roughly  52.  This  coarse  mesh  accuracy  more  than  makes  up  for  the  factor  of  3 
difference  in  computational  cost  between  the  EC  and  LC  quadratures. 

Lastly,  we  remark  on  some  issues  encountered  in  producing  our  coarse 
mesh  results.  In  initial  attempts  to  solve  this  problem,  the  EC  method  failed 
to  converge.  This  type  of  error  had  been  encountered  before  in  exploring  set- 
to-zero  fix-ups  for  the  source  moments.  However,  there  are  no  set-to-zero  fix¬ 
ups  in  TETRAN.  We  examined  the  root-solver  and  found  that  the  error  was 
occurring  because  of  an  inconsistent  implementation  of  the  transition  from 
the  asymptotic  solution  to  the  source  system  (requiring  no  root-solves)  to  the 
first  guess  algorithm  and  the  root-solver  (discussed  in  Appendix  B).  This 
caused  the  source  to  be  calculated  inconsistently  between  iterations  and 
prevented  convergence.  The  problem  was  fixed  by  changing  the  break  points 
to  be  moderately  more  conservative  regarding  the  use  of  the  asymptotic 
formulas.  Thus,  the  root-solver  is  invoked  more  often  in  these  extreme  cases, 
producing  consistently  accurate  source  moments.  At  some  point  in  the  future, 
we  will  have  to  return  to  this  issue  and  re-evaluate  the  region  of  applicability 
of  the  asymptotic  solutions.  However,  the  performance  penalty  for  the  more 
conservative  approach  is  minimal  and  thus  does  not  warrant  a  high  priority 
for  future  work. 

Multi-group,  Anisotropic  Scattering  Performance 

Raising  the  level  of  complexity  of  our  test  problems  higher,  we  examine 
TETRAN's  ability  to  solve  multi-group  problems  with  and  without  anisotropic 
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scattering.  The  problem  examined  is  similar  to  one  found  in  the  literature 
(Castrianni,  1997).  We  call  it  the  water-iron-water  problem. 

The  water-iron- water  problem  is  comprised  of  three  regions.  Region  1 
is  a  20  centimeter  cube  of  water  centered  at  the  origin  with  a  uniformly 
distributed  source  of  particles  emitted  isotropically.  Region  2  is  a  hollow  iron 
cube  that  surrounds  the  source  and  has  a  thickness  of  10  centimeters. 
Finally,  region  3  is  a  hollow  cube  filled  with  water  that  surrounds  the  iron 
region  and  has  a  thickness  of  30  centimeters.  This  geometry  is  depicted  in 
Figure  21. 


Figure  21.  Geometry  for  the  Water-Iron-Water  Problem. 

This  problem  is  driven  by  a  mono-energetic  source  of  strength  1.0 
particles/sec  emitted  in  group  one.  Vacuum  boundary  conditions  are  assumed. 
Two  anisotropic  scattering  approximations  were  used.  For  the  EC  method,  we 
used  Pi  scattering  cross  sections  that  were  corrected  using  the  diagonal 
transport  approximation  (MacFarlane,  1992).  Thus,  we  get  Pi  accuracy  from 
an  isotropic  cross  section  calculation.  These  cross  sections  were  also  used  for 
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an  LC  calculation.  Additionally,  we  solved  the  problem  using  LC  and  P3 
scattering  moments  that  were  corrected  using  the  diagonal  transport 
approximation  to  investigate  the  effect  and  performance  of  our  anisotropic 
scattering  treatment.  We  attempted  to  solve  the  P3  problem  using  EC  but 
were  unable  to  do  so  because  using  Legendre  moments  and  spherical 
harmonics  to  update  the  scattering  source  produces  negative  sources.  EC 
requires  positive  source  data  in  order  that  we  may  calculate  the  source 
coefficients.  This  is  especially  true  of  water  where  the  scattering  off  the 
hydrogen  atom  is  particularly  forward  peaked  because  neutrons  are  only 
forward  scatter  in  the  laboratory  frame  of  reference.  Attempting  to  use 
Legendre  polynomials  to  fit  such  a  function  invariably  produces  negative 
values.  Additionally,  narrow  energy  groups  can  cause  negative  cross  sections 
to  be  produced  for  certain  materials.  Alternative  cross  section  representations 
will  be  needed  in  order  to  eliminate  negative  cross  sections.  This  is  a  subject 
for  further  research  (DelGrande,  1998).  The  first  three  neutron  energy  groups 
(12.0  MeV  to  17.0  MeV)  in  the  MATXS10  cross  section  library  (MacFarlane, 
1992)  were  used  for  this  problem.  They  are  presented  in  Table  14,  Table  15, 
and  Table  16. 

The  MATXS10  cross  sections  were  selected  for  this  research  for  two 
reasons.  They  are  readily  available  in  the  TRANSX  package  and  contain 
many  common  engineering  materials;  they  include  up  though  P4  scattering 
moments;  and  they  are  based  on  the  ENDF/B-VT  evaluations  (the  most 
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current  nuclear  data  files).  Also,  MATXSlO's  energy  structure  is  identical  to 
that  of  the  MCNP  multi-group  cross  section  library,  which  is  based  upon  the 
MENDF5  library  used  at  Los  Alamos  for  high  energy  transport  calculations. 
However,  there  are  substantial  differences  between  the  MATXS10  cross 
sections  and  their  MCNP  equivalent.  First,  the  MATXS10  library  is  based  on 
ENDF/B-VI  while  the  MCNP  version  is  derived  from  ENDF/B-V  data.  More 
importantly,  in  order  to  be  incorporated  into  MCNP,  the  multi- group  P0 
through  P4  data  were  processed  into  equally-probable  cosine  bins  to  treat  the 
scattering  using  the  MCNP  data  structure  (MCNPXS,  1997).  The  data  must 
also  have  been  processed  to  eliminate  the  possibility  of  negative  cross 
sections,  which  is  as  intolerable  for  the  Monte  Carlo  method  as  it  is  to  EC. 
Thus,  although  the  cross  sections  appear  to  be  derived  from  a  common 
parent,  they  are  at  least  partly  different  because  of  the  data  requirements  of 
the  codes  for  which  they  are  intended.  Because  of  the  substantial  differences 
between  multi-group  and  continuous  energy  cross  sections,  we  compare 
TETRAN's  results  with  MCNP's  multi-group  results.  This  at  least  minimizes 
the  differences  due  to  the  cross  section  representations.  Of  course,  the  MCNP 
results  do  not  have  the  angular  and  spatial  discretization  errors  that  are 
inherent  in  the  TETRAN  results. 
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Table  14.  Pi  Transport  Corrected  Cross  Sections  for  Water. 


Water 

cfs>g'_>g(cm  x) 

g 

CTg(cm  x) 

g'  =  1 

g'  =  2 

g'  =  3 

1 

0.074512 

2.07399E-03 

2 

0.079409 

1.31405E-02 

1.75144E-03 

3 

0.081262 

5.74009E-03 

1.42416E-02 

2.85886E-03 

Natural 

Iron 

oS)g^g(cm'1) 

g 

ag(cm_1) 

g'  =  l 

g'  =  2 

g'  =  3 

1 

0.132931 

1.08255E-02 

2 

0.137572 

1.38718E-02 

1.26268E-02 

3 

0.138529 

1.41291E-03 

1.34612E-02 

1.74764E-02 

Table  15.  P3  Transport  Corrected  Cross  Sections  for  Water. 


1=0 

CTs/,g'->g(cm 

g 

og(cm  *) 

g'  =  l 

g'  =  2 

g'  =  3 

1 

0.084473 

1.20352E-02 

2 

0.089661 

1.31405E-02 

1.20036E-02 

3 

0.092259 

5.74009E-03 

1.42416E-02 

1.38556E-02 

Z=1 

CTsZ,g'^g(cm_1) 

g 

g'  =  1 

g'  =  2 

g'  =  3 

1 

9.96124E-03 

2 

8.13291E-03 

1.02522E-02 

3 

2.69041E-03 

7.40822E-03 

1.09968E-02 

1=2 

<TsZ,g'^g(cm~1) 

g 

g'=l 

g'  =  2 

g'  =  3 

1 

6.67593E-03 

2 

3.61894E-03 

7.21434E-03 

3 

2.54424E-03 

2.18526E-03 

6.84057E-03 

1=3 

°s!,g’->g(cni"1) 

g 

g'=l 

g'  =  2 

g'  =  3 

1 

3.18798E-03 

2 

2.60006E-03 

3.58607E-03 

3 

2.01574E-03 

2.35980E-03 

3.02045E-03 
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Table  16.  P3  Transport  Corrected  Cross  Sections  for  Natural  Iron. 


1=0 

as/,g'-»g(cm  x) 

g 

cg(cm  *) 

g'  =  1 

g'  =  2 

g'  =  3 

1 

0.162578 

4.04721E-02 

2 

0.173704 

1.38718E-02 

4.87585E-02 

3 

0.180036 

1.41291E-03 

1.34612E-02 

5.89831E-02 

1=1 

os;,g^g(cm'1) 

g 

g'  =  l 

g'  =  2 

g'  =3 

1 

2.96467E-02 

2 

3.08452E-03 

3.61317E-02 

3 

1.58912E-04 

9.05030E-04 

4.15068E-02 

1=2 

as/ig^g(cm_1) 

g 

g'  =  l 

g'  =  2 

g'  =  3 

1 

1.81862E-02 

2 

1.00274E-03 

2.24502E-02 

3 

1.79372E-04 

1.03198E-03 

2.75191E-02 

1=3 

as/,g'->g(cm_1) 

g 

g'  =  l 

g'  =  2 

g'  =  3 

1 

8.75387E-03 

2 

7.57219E-04 

1.12032E-02 

3 

9.07334E-07 

4.47590E-04 

1.39011E-02 

The  results  of  TETRAN's  multi-group  transport  calculations  are 
presented  below.  The  data  presented  are  for  three  different  tetrahedral 
meshes:  528,  8849,  and  17369  cells.  Table  17  presents  the  group  dependent 
average  cell  optical  thickness  by  region.  This  is  followed  by  the  results 
presented  by  region. 
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Table  17.  Group  and  Mesh  Dependent  Cell  Optical  Properties  for  Problem  6. 


Region 

Cells 

gEi 

eg 

g=2 

gE 3 

Source: 

48 

0.64718 

0.68971 

0.70581 

162 

0.41489 

0.44216 

0.45248 

161 

0.42145 

0.44915 

0.45963 

Iron: 

144 

1.4830 

1.5348 

1.5454 

931 

0.78105 

0.80832 

0.81394 

1350 

0.68926 

0.71332 

0.71828 

Water: 

336 

1.6522 

1.7608 

1.8019 

7756 

0.57544 

0.61326 

0.62757 

15858 

0.44951 

0.47906 

0.49023 

LC  EC  LC 

(P,)  <p,)  (Pa) 


528  cells 

i _ ; _ 1 

oo4i7  COlls 

17369  cells 

.  .  . 

MCNP  High 

MCNP  Mean 

.  .  . 

MCNP  Low 

Figure  22.  Group  1  Scalar  Flux  in  Source  Region. 
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Figure  23.  Group  2  Scalar  Flux  in  Source  Region. 
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The  above  figures  show  that  there  is  a  substantial  difference  in  results 
between  using  the  Pi  and  P3  scattering  approximations.  They  also  show  good 
agreement  between  the  LC  P3  and  the  MCNP  P4  calculation  (indeed,  the  LC 
P3  cross  sections  are  modified  P4  cross  sections).  The  group  1  flux  shows 
excellent  agreement  with  MCNP.  However,  we  should  expect  this  because 
group  1  is  the  highest  group  for  this  problem.  It  is  in  the  downscatter  groups 
that  we  expect  to  see  the  effect  of  the  different  treatments  of  anisotropic 
scattering.  Indeed,  we  see  that  both  group  2  and  group  3  scalar  fluxes  are 
substantially  different  from  MCNP  for  the  Pi  scattering  cross  sections 
whereas  the  LC  P3  calculation  is  close  to  the  MCNP  results  for  all  three 
groups.  We  now  present  the  iron  region  results. 


Figure  25.  Group  1  Scalar  Flux  in  Iron  Region. 
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The  discrete  ordinates  results  for  the  scalar  flux  in  the  iron  region 
show  excellent  agreement  with  the  MCNP  results  for  both  the  Pi  and  P3 
calculations.  This  is  because  neutron  scattering  in  iron  is  more  isotropic  than 
in  water  due  to  its  greater  nuclear  mass.  Another  interesting  aspect  of  the 
iron  region  is  that  the  EC  calculation  is  more  accurate  than  the  LC  solution 
as  compared  to  MCNP  for  all  three  energy  groups.  The  coarse  mesh  LC 
calculation  is  measurably  off  for  the  fluxes  but  converges  for  the  finer 
meshes.  This  is  because  the  average  optical  thickness  of  a  cell  in  the  iron 
region  is  1.52  (over  all  groups)  for  the  coarse  mesh  and  only  0.80  and  0.71  for 
the  two  finer  meshes,  respectively.  Additionally,  the  maximum  cell  optical 
thicknesses  range  from  2.7  for  the  coarse  mesh  to  1.7  for  the  fine  mesh  (for  all 
energy  groups).  Thus,  we  see  again  the  thick  cell  performance  of  EC  as 
compared  to  LC.  We  now  will  examine  the  scalar  flux  results  in  the  outer 
water  region. 
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Figure  28.  Group  1  Scalar  Flux  in  Outer  Water  Region. 
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Figure  29.  Group  2  Scalar  Flux  in  Outer  Water  Region. 
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Figure  30.  Group  3  Scalar  Flux  in  Outer  Water  Region. 

We  see  the  same  type  of  behavior  for  the  water  region  that  was  seen 
for  the  source  region.  Both  Pi  calculations  are  substantially  different  from  the 
LC  P3  calculation.  However,  the  LC  P3  calculation  compares  very  well  to  the 
MCNP  calculation.  However,  this  should  be  expected  because  water  is  a 
pronouncedly  anisotropic  scatterer  of  neutrons.  Notice  also  the  coarse  mesh 
performance  of  EC  compared  to  LC.  In  this  instance,  the  average  cell  optical 
thickness  (over  all  groups)  is  1.74  for  the  coarse  mesh  and  only  0.61  and  0.47 
for  the  finer  meshes.  The  maximum  optical  thicknesses  range  from  3.5  for  the 
coarse  mesh  to  1.0  for  the  fine  mesh  (nearly  constant  over  all  groups).  Notice 
also  that  the  group  1  fluxes  compare  favorably  with  the  MCNP  results  and 
the  variations  occur  in  the  lower  groups  as  expected. 
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Finally,  we  examine  the  group  currents  emerging  from  the  +X  face  at 
X=50.0  cm.  The  behavior  of  the  current  at  the  boundary  for  this  problem 
gives  a  measure  of  the  impact  of  the  various  methods  and  Legendre 
scattering  order  (Pi  or  P3)  on  the  particle  leakage. 


Figure  31.  Group  1  Current  Out  of  +X  Face  (X=50.0  cm). 
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Figure  32.  Group  2  Current  Out  of  +X  Face  (X=50.0  cm). 
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Figure  33.  Group  3  Current  Out  of  +X  Face  (X=50.0  cm). 
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From  the  above  figures,  it  is  obvious  that  the  current  is  particularly 
sensitive  to  the  scattering  approximation.  The  Pi  calculations  are  not  very 
good  for  any  group  in  this  case.  The  P3  calculation  is  much  better.  In  order  to 
make  sure  that  the  algorithm  was  working  correctly,  particle  balance 
calculations  were  done,  which  showed  that  particles  are  being  conserved. 
This  forces  us  to  conclude  that  the  differences  between  the  current 
calculation  are  due  to  the  nuclear  data  and  the  methods  (Monte  Carlo  vs. 
discrete  ordinates  vs.  EC/LC  approximation)  used.  This  conclusion  could  be 
further  validated  if  a  Monte  Carlo  code  existed  which  forced  particles  to 
travel  along  discrete  ordinates  angular  quadrature  directions  and  was 
multigroup  capable.  Such  a  code  would  produce  results  which  contained 
angular  and  energy  discretization  errors  but  would  not  have  spatial 
dicretization  error,  making  it  a  superior  benchmarking  code  for  discrete 
ordinates  results.  The  initial  results  of  such  a  code,  MCSN,  were  presented 
during  testing  of  the  convergence  of  EC  and  LC  in  problem  one.  Currently,  no 
such  capability  exists  for  multigroup,  anisotropic  scattering  problems.  Dr. 
Kirk  Mathews  (AFIT)  is  developing  an  upgrade  to  MCSN  to  make  it 
multigroup  capable  and  able  to  use  the  T-matrix  data  structure  (Mathews, 
1998).  It  will  be  interesting  to  see  how  EC  performs  once  positive  T-matrix 
data  becomes  available  for  this  problem  (DelGrande,  1998). 

The  above  figures  show  TETRAN  is  performing  multi-group, 
anisotropic  scattering  calculations  accurately  (given  the  correct  data).  The 
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optical  thickness  of  the  regions  between  the  8849  cell  mesh  and  the  17369  cell 
meshes  are  not  very  different.  This  is  true  of  the  source  region  as  well.  These 
trends  are  indicative  of  other  problems  that  current  mesh  generators  have 
with  regard  to  mesh  refinement.  The  mesh  generator  used  in  this  work 
produced  meshes  that  were  inconsistent  with  increasing  refinement,  i.e.  finer 
meshes  were  produced,  but  not  all  of  the  regions  in  the  problem  were  refined 
consistently.  This  behavior  along  with  the  lack  of  volume  conservation  raises 
concerns  regarding  the  use  of  common  CAD/CAM  codes  for  transport 
calculations.  The  T-matrix  implementation  for  the  group  to  group  scattering 
source  appears  to  function  well,  as  does  the  iteration  approach  for  each 
energy  group.  LC  is  theoretically  capable  of  handling  any  order  of  anisotropic 
scattering  while  EC  will  need  positive  cross  sections  before  it  will  be  usable 
for  problems  requiring  anisotropic  scattering. 

HPF  Parallel  Demonstration 

In  this  last  test,  we  demonstrate  the  parallel  scaling  performance  and 
accuracy  of  our  HPF  parallel  implementation  of  TETRAN.  For  this  problem, 
we  use  the  162  cell  mesh  used  in  the  convergence  rate  problem  presented 
earlier.  The  results  are  for  both  EC  and  LC.  We  chose  to  use  this  simple 
problem  because  the  HPF  compiler,  PGHPF  2.4  (development  compiler), 
produced  code  that  runs  between  6  (EC)  and  8  (LC)  times  slower  than  the 
IBM  xlf90  serial  code  used  in  the  tests  presented  earlier.  We  chose  a 
relatively  small  mesh  because  ASC's  IBM  SP  proved  to  have  problems  with 
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its  message  passing  library  that  prevented  us  from  using  more  than  two 
processors  for  any  larger  problem. 

For  our  initial  evaluation  of  TETRAN's  parallel  performance,  we  look 
at  the  consistency  of  the  answers  and  the  scaling  of  the  code  with  increasing 
processors.  We  look  at  consistency  to  ensure  that  the  HPF  code  produces  the 
same  results  regardless  of  the  number  of  processors  used  (deterministic).  The 
speed  up,  S,  is  used  to  determine  the  parallel  scaling  of  the  code.  Speed  up  is 
the  ratio  of  the  time  to  run  the  problem  on  one  processor  to  that  for  multiple 
processors.  No  further  investigation  was  done  because  the  compiler  has  too 
many  performance  issues  to  make  any  further  discussions  useful.  Future 
investigations  should  look  at  compiling  with  more  than  one  HPF  compiler  if 
available. 

The  HPF  code  was  found  to  provide  the  same  results  regardless  of  the 
number  of  processors  used  as  required  by  the  HPF  standard  (Koelbel,  1994). 
For  LC,  the  scalar  flux  was  1.541295  particles/(cm3-s).  The  EC  calculation  for 
the  scalar  flux  was  1.540445  particles/(cm3-s).  The  parallel  speed-up  (S)  and 
efficiency  (E)  are  presented  in  Table  18. 
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Table  18.  HPF  Parallel  TETRAN  Results. 


Run  Time 
(s) 

Speed-up,  S 

Efficiency,  E 
(S/N) 

Nproc 

LC 

EC 

LC 

EC 

LC 

EC 

l 

882.93 

2833.8 

1 

1 

1.000 

1.000 

2 

451.35 

1435.6 

1.9562 

1.974 

0.978 

0.987 

4 

256.43 

752.87 

3.4431 

3.764 

0.861 

0.941 

8 

154.73 

424.31 

5.7063 

6.6786 

.  0.713 

0.835 

10 

114.92 

321.99 

7.6831 

8.8009 

0.768 

0.880 

16 

110.96 

327.43 

7.9571 

8.6547 

0.497 

0.541 

32 

170.32 

248.33 

5.1838 

11.412 

0.162 

0.357 

1  (xlfOO) 

107.74 

489.31 

N/A 

N/A 

N/A 

N/A 

We  see  in  Table  18  that  the  HPF  code  is  scaling  to  about  8  processors  for  LC 
and  about  10  processors  for  EC  before  saturating.  Beyond  16  processors,  both 
methods  appear  to  saturate  as  evidenced  by  the  drop  in  the  parallel  efficiency 
shown  in  Table  18.  Using  32  processors  is  markedly  worse  for  LC  than  using 
8  processors.  The  EC  method  performs  better  than  the  LC  method, 
presumably  because  the  processors  are  kept  busy  performing  the  expensive 
quadrature  as  was  hoped.  The  saturation  is  probably  due  to  the 
interprocessor  communications  required  in  the  source  and  scalar  flux  updates 
or  a  demonstration  of  Amdahl's  Law  which  basically  says  that  the  part  of  a 
program  that  can't  be  optimized  (in  our  case,  the  source  updates)  will 
eventually  dominate  the  runtime  for  a  code  (Dowd,  1998).  Because  the 
compiler  has  performance  issues,  we  can't  be  sure.  Clearly,  much  further 
investigation  is  required.  A  plot  of  the  speed  up  versus  the  number  of 
processors  is  shown  in  Figure  34. 
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.  Number  of  Processors 


Figure  34.  Demonstration  of  Speed  Up  for  HPF  Parallel  TETRAN. 

Regardless  of  the  above  issues,  it  is  encouraging  that  the  code  compiles 
and  runs  and  that  the  EC  method  appears  to  scale  better  than  LC.  This 
validates  our  design  methodology  of  developing  TETRAN  using  standard 
Fortran  90  and  then  inserting  HPF  directives  as  needed.  We  worked  very 
closely  with  The  Portland  Group  (Larry  Meadows)  to  work  out  the  bugs  that 
prevented  the  compilation  of  TETRAN.  Portland  Group  is  now  working  on 
further  enhancing  their  compiler  to  get  the  performance  that  is  readily 
available  with  the  serial  Fortran  90  compiler  available  on  the  IBM  SP.  They 
believe  that  the  problems  with  the  compiler  are  in  their  implementation  of 
the  Fortran  90  intrinsic,  MATMUL.  If  this  is  the  case,  then  the  performance 
results  we  report  are  preliminary  because  we  use  the  MATMUL  intrinsic  in 
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every  computational  subroutine  in  TETRAN.  It  is  also  at  the  core  of  the 
source  update  routines.  We  expected  that  the  maturity  of  HPF  compilers 
would  be  better  than  was  found  (PGHPF  was  supposed  to  be  mature). 
However,  since  HPF  is  still  a  developing  compiler  area,  providing  transport- 
related  challenges  to  the  compiler  developers  will  hopefully  results  in  better 
HPF  compilers  for  the  future.  Lastly,  in  addition  to  the  compiler  issues,  we 
found  that  the  IBM  SP  has  problems  running  the  HPF  version  of  TETRAN 
for  large  problems.  This  problem  manifests  itself  as  message-passing  errors. 
The  staff  at  the  ASC  MSRC  is  working  on  this  issue  as  other  users  have 
complained  of  it  as  well.  Although  future  efforts  should  consider  using  the 
IBM  SP  because  of  its  capability  to  handle  large  amounts  of  data,  other 
capable  machines  are  available  at  the  MSRC.  Most  noteworthy  is  the  SGI 
ORIGIN  2000.  This  machine  is  ASC's  largest  shared  memory  machine  (ASC, 
1998).  It  is  possible  that  the  HPF  version  of  TETRAN  will  perform  better  on 
such  a  machine  versus  the  IBM  SP.  This  performance  difference  has  been 
observed  by  other  researchers  using  different  codes  (Little,  1998). 

This  chapter  presented  TETRAN  test  results.  Each  problem  tested 
critical  features  of  the  spatial  quadratures  and  the  overall  implementation  of 
the  code.  The  results  show  that  TETRAN  is  performing  as  expected.  Both 
spatial  quadratures  converge  to  the  same  solution  for  the  same  problem  and 
appear  to  be  accurate  compared  to  MCNP  calculations.  The  EC  quadrature 
proved  to  be  extremely  robust  as  evidenced  by  the  problems  it  solved.  The  T- 
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matrix  approach  accurately  solved  the  multigroup,  anisotropic  scatter 
problem  and  aided  in  the  parallelization  of  TETRAN.  Three  major  issues 
(beyond  the  scope  of  this  research)  were  uncovered  during  this  testing 
regarding  mesh  generation,  positive  nuclear  data,  and  HPF  compilers. 
However,  TETRAN  proved  to  be  a  strong  performer  on  a  variety  of  difficult 
problems. 
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Chapter  V:  Conclusions  and  Recommendations 
Conclusions 

We  derived  a  new  form  of  the  EC  spatial  quadrature  that  is  both 
elegant  and  illuminating  with  regard  to  its  numerical  implementation.  This 
new  form  of  the  EC  quadrature  relies  on  the  use  of  direct  affine 
transformations  to  pass  cell  and  flux  moments  within  a  mesh.  EC  and  LC 
were  implemented  into  a  transport  code,  TETRAN.  Both  methods  are 
numerically  robust,  convergent,  and  accurate.  Several  numerically  accurate 
and  robust  algorithms  (for  MQ{x,y,z,w)  and  the  ^-functions)  were  developed, 
enabling  accurate  computation  of  the  EC  (Appendix  A)  and  LC  quadratures. 
We  tested  the  method's  performance  on  a  variety  of  problems,  demonstrating 
convergence  and  robustness  for  poorly  shaped  and  optically  thick 
computational  cells.  The  EC  spatial  quadrature  is  computationally  expensive, 
costing  approximately  3  ms/phase  space  cell  while  LC  costs  about  a  third  of 
‘that  at  1  ms/phase  space  cell. 

A  new  root-solver  was  developed  that  uses  a  very  accurate  first  guess 
algorithm  based  upon  the  characteristics  of  the  solution  space  (Appendix  B). 
The  capability  to  perform  standard  multi-group  calculations  with  anisotropic 
scattering  was  demonstrated  for  a  difficult  problem  with  great  success.  Our 
multi-group  implementation  uses  the  T-matrix  approach  over  the  traditional 
angular  flux  moments  ( )  approach  for  updating  the  scattering  source. 
However,  use  of  EC  on  greater  than  first  order  anisotropic  scattering 
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problems  will  have  to  wait  until  research  is  complete  on  a  positive  cross 
section  generator.  Standard  cross  sections,  using  spherical  harmonics  with 
Legendre  moments  of  the  scattering  cross  section,  produce  negative  sources 
which  are  intolerable  (and  physically  wrong)  for  EC.  The  LC  method  did  not 
have  any  problems  using  the  standard  data  since  the  LC  source  assumptions 
do  not  require  positivity  of  the  source  moments.  Finally,  we  demonstrated  the 
parallel  operation  of  an  HPF  version  of  TETRAN  on  a  simple  problem.  Both 
EC  and  LC  scale  with  increasing  processors,  the  EC  method  scaling  better 
than  LC.  Future  work  in  this  area  will  have  to  wait  for  a  more  mature 
compiler. 

In  summary,  this  effort  developed  a  parallel,  unstructured  tetrahedral 
mesh  radiation  transport  code,  TETRAN.  This  code  solves  the  BTE  using 
either  the  LC  or  the  EC  method  using  direct  affine  transformations  to 
calculate  the  needed  source  coefficients.  It  can  solve  multigroup,  anistropic 
scattering  problems  using  a  new  approach  to  calculating  the  scattering  source 
—  the  T-matrix  algorithm.  New  findings  regarding  the  nature  of  the  source 
and  inflow  flux  equations  allowed  the  implementation  a  new  non-linear  root¬ 
solving  strategy.  HPF  directives  made  the  code  parallel  on  the  IBM  SP. 
Substantial  effort  was  made  to  implement  the  EC  method  such  that  it  was 
robust  and  accurate.  Seven  test  problems  demonstrate  this  robustness  and 
accuracy  as  discussed  earlier. 
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Recommendations  for  Future  Work 

Several  potential  areas  of  future  research  were  uncovered  during  this 
research:  mesh  generation,  parallel  algorithms,  and  benchmarking. 

If  unstructured  tetrahedral  mesh  transport  codes  are  to  succeed,  they 
must  have  good  meshes  on  which  to  operate.  We  recommend  that  a  new  mesh 
generator,  optimized  for  radiation  transport,  be  developed.  It  should 

-  conserve  volumes, 

-  minimize  aspect  ratios  of  cells, 

-  provide  users  control  of  optical  thicknesses, 

-  combine  triangular  cells  on  selected  surfaces  with  tetrahedral  cells  that 
align  with  them,  and 

-  provide  optional  self-similar  mesh  refinement  where  feasible. 

We  recommend  that  MCSN  be  extended  to  multigroup,  anisotropic 

scatter  using  T-matrices,  in  order  to  provide  more  useful  benchmark  results. 

Current  computer  architectures  re  dependent  upon  locality  of  reference 
for  effective  use  of  caches.  We  recommend  development  of  a  regular- 
hexahedral-mesh  code  analogous  to  TETRAN. 

Various  enhancements  are  needed  to  make  TETRAN  more  useful, 
including: 

-  generating  positive  T-matrices,  and 

-  adding  convergence  acceleration. 

Other  enhancements  might  improve  the  speed  of  the  code,  such  as: 
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-  using  balance  equations  to  eliminate  expensive  moment  functions  where 
well-conditioned,  and 

-  exploring  alternatives  for  indirect  addressing  of  fluxes. 

HPF  and  Fortran  90/95  compilers  are  not  yet  mature.  Many  are  still  in 

early  stages  of  development.  We  recommend  testing  different  hardware 
platforms  and  compilers.  Ultimately,  it  may  be  necessary  to  recode  using  MPI 
and/or  using  various  coding  tricks  to  squeeze  performance  from  a  specific 
platform  and  compiler.  We  recommend  that  this  be  postponed  until 
everything  else  is  done,  both  to  avoid  wasted  effort  and  because  it  may  prove 
to  be  unnecessary. 
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Appendix  A:  Exponential  Moment  Functions  —  Evaluation  and 

Identities 

The  EC  method  requires  the  evaluation  of  a  number  of  exponential 
moment  functions.  This  appendix  presents  several  moment  function 
identities  used  in  the  derivation  of  EC.  The  algorithm  for  numerically 
evaluating  four  argument  exponential  moment  functions  is  outlined. 

Definitions 

As  previously  presented,  the  exponential  moment  functions  are  defined 
as 

^ ni—1 

Mn(x1,x2,...,xm)  =  jdtj(l- tj)n  e~Xl ll | dt2  e(xi~Xs)t2 . . .  J dtm  e(Xm-t~Xm)t'n .  (82) 

o  oo 

Note  that  A1n(x1,x2,...xm)  >  0.  Additionally,  the  moments  functions  are 
orderless  with  respect  to  their  arguments,  e.g. 

Atn(x1,x2,...xm)  =  Afn(x2,x1,...xm)  and  so  on.  The  multi-argument  Mq 
function  is  most  often  used  in  the  EC  method.  The  ratio 
Ali(x1,x2,...,xm)/  At0(x1,x2,...,xm)  also  appears  frequently  with  respect  to 
EC.  Thus,  we  define 


p(x1,x2,...xm)  = 

In  this  case,  0  <  p(x1,x2,...xm)  <1, 
A^and  Mq. 


M1(x1,x2,..,xm) 

M0(x1,x2,...xm)' 

which  follows  from  the  definition  of 


(83) 


117 


Lastly,  we  define  a  new  exponential  moment  function  ratio,  *Rj ,  given 


as 


7?j(xi,x2 


Af0(x1,x2,...,xm,xj) 

M)(x1,x2,...,xm) 


(84) 


where  the  j  subscript  indicates  that  the  jth  argument  (pi,  2,  ....or  m)  is 
repeated  in  the  numerator.  Note  that  the  numerator  of  *Rj  is  a  function  of 

m+1  coefficients  whereas  the  denominator  is  a  function  of  m  terms  and  that 
0<^j(x1x2,...,xm)<l. 


Exponential  Moment  Function  Identities 

Many  exponential  moment  function  identities  are  used  in  the  EC 
spatial  quadrature  derivation.  These  identities  are  used  to  derive  the 
quadrature  and  cast  it  in  a  stable  form.  Most  were  presented  in  (Mathews, 
1997).  All  of  the  presented  identities  were  verified  symbolically  using 
Mathematica®  3.0.  Additional  detail  regarding  exponential  moment 
functions  can  be  found  in  (Minor,  1993). 


Identity  1 


A4n(x1,x2,...,xm)  = 


_  '^4n(xj,...,Xj_^,Xj+2,...,xm)  A4n(x2,...,xj£_1,xk+1,...,xm) 


Xj-Xk 


(85) 


Xj  *Xk 
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The  divided  difference  identity  is  perhaps  the  oldest  of  our  identities.  It  was 
initially  proven  in  Minor's  dissertation  (Minor,  1993).  The  above  identity  is 
used  to  stably  evaluate  ^j(x1,x2,...,xm)  and  p(x1,x2,...,xm) . 


Identity  2 

Mo(x1,x2,...,xm)  =  exp(-x1)M0(-x1,x2  -x1(...,xm  -xx)  (86) 


Identity  2  is  used  to  cast  <7?j(x1,x2,...,xm)  in  a  numerically  stable  form.  It  is 
valid  for  all  values  of  m. 


Identity  3 


Mn(Xi,...,Xk 


>xm )  ■^dn(xj,...,xk,xk,...,xm) 


(87) 


Identity  3  is  used  in  the  evaluation  of  the  source  system  and  inflow  face  flux 
system  Jacobians  as  well  as  in  the  derivation  of  the  EC  spatial  quadrature.  It 
was  verified  for  all  n  and  m. 


Identity  4 

Ai0(0,x1,x2,...,xm)  =  Af1(x1,x2,...,xm)  (88) 

Identity  4  is  used  with  Identity  2  to  recast  <Rj(x1,x2,...xm)  into  a  numerically 
stable  form  involving  p(xx  ,x2 ,. . . ,xm ) .  The  identity  was  verified  for  all  values 
of  m. 
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Identity  5 


M0(x1,x2,...,xm) 


Af0(x1,x2,...,xm_1)-exp(-x1)A10(x2  -x1,...,xm  -x^ 

xm  (89) 

xm 


Identity  5  is  used  in  the  numerical  evaluation  of  the  moment  functions.  It 
was  confirmed  by  direct  symbolic  evaluation  of  the  integral  functions  using 
Mathematica®  for  m=2  to  5.  This  identity  was  first  introduced  in  (Mathews, 
1996). 


Identity  6 

m 

p(x1,x2,...,xm)  +  ]Tftj(x1,x2)...,xm)  =  1  /  (90) 

j=i 

Identity  6  is  new  to  this  research.  It  is  used  to  recast  the  EC  spatial 
quadrature  into  a  numerically  stable  form  without  subtractions.  It  was 
confirmed  by  direct  symbolic  evaluation  of  the  integral  functions  using 
Mathematica®  for  m=2  to  4. 

A1o(x,y,z,w)  and  p(x,y,z,w)  Numerical  Evaluation 

The  development  of  stable  and  accurate  routines  to  evaluate  the 
exponential  moment  functions  and  their  ratios  was  an  important  contribution 
of  this  research.  Previously,  Mathews  developed  streamlined  routines  to 
evaluate  Afo(x)  and  p(x),  A1o(x,y)  and  p(x,y),  and  AIo(x,y,z) and  p(x,y,z). 
The  algorithms  are  presented  in  (Mathews,  1997).  For  this  effort,  we 
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extended  the  algorithm  used  for  three-argument  case  to  the  four-argument 

case.  Although  the  four-argument  algorithm  is  not  as  elegant  as  the  three- 

argument  one,  it  is  stable  and  accurate  to  at  least  10  digits.  The  pseudo-code 

for  the  four-argument  function  algorithm  is  presented  below. 

Sort  x,  y,  z,  and  w  into  ascending  order  -  w>z>y>x.  We  can  do  this 
because  the  exponential  moment  functions  are  orderless  with  respect  to 
their  arguments. 

If  x>38.0  then  all  arguments  are  large  enough  that 

kA  (  ,1  ,  N  1111 

M0(x/y/z/w)  = - and  p(x,y,z,w)  =  l - .  This 

xyzw  x  y  z  w 

approximation  is  accurate  to  at  least  12  digits. 

If  abs(x)  and  abs (w)  <  0.025,  (all  arguments  small)  then  we  use  a 
Horner  nested  MacLaurin  series  where  we  expand  both  A1o(x,y  #z  ,w)  and 
p(x,y#Z,w)to  5ch  order  terms  in  x,  y,  z,  and  w  and  retain  all  terms  up 


to 

5th  order. 

y-x 

z 

~x 

Let 

x  = 

and  (0  : 

=  — 

- . 

w-x 

w 

-x 

If 

X  <  0.1  or 

X  >  0.9 

or 

co  <  0.1  or  CO  >  0.9  or  (co  —  x)  <  0.1  or  (co  —  x)  >  0.9 

or  abs(w-x)<  0.1  Max{abs  (x)  ,  abs  (w)  }  then  one  or  more  arguments  are 
close  together  but  not  near  zero.  In  this  case,  we  use  identity  5 
(above)  to  get  : 

M0(XfYiZiW)=  -exp(-x)A^(y-x,z-x,w-xj 

w 

(x  v  =  Mo(x>y>*)p(x>y>*)-Mo(x,y,z,w) 
w  Mo(x,y>z,w) 

(Note  that  if  abs (x)  >  abs (w) ,  we  switch  x  and  w  for  better 
conditioning. ) 

Else,  we  use  the  general  divided  difference  formulation  given  as: 

Mq ( X  ,y ,  z ,  w)  =  [  (1  -  x)((0  -  x)(l  -  ©  )Mq  (x)  -  (1  -  CO )  CO  Ato  (y)  + 
x(i-x)M)(z)  -x®(®-x)M>(w)]/ 

[{w-xf  (1-X)X(®~X)(1-®)®] 
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p(x/y,z,w)=[(l-x)(©-x)(l-©)p(x)  A1q  (x)  -  (1  -  co  )  co  p(y )  Mq  (y )  + 
(1-X)XP(Z)M)(Z)  -X©(©-X)  p(w)Mo(w)]/ 

[(1  ~  X)  (©  -  X)  (1  -  © )  Mq  (x)  -  (1  -  to )  ©  Mq  (y)  + 

(i-x)x  Mb(z)-X®(©~X)  M)(w)] 

The  last  equations  are  obtained  by  applying  Identity  1  to  A\q  and  Mi , 
rearranging  terms,  and  applying  the  definitions  for  %  and  to  . 
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Appendix  B:  Evaluation  of  Source  and  Inflow  Flux  Coefficients 
The  exponential  characteristic  method  requires  that  source  and  inflow 
flux  coefficients  be  computed  and  used  in  the  sub-cell  spatial  quadrature.  In 
this  appendix,  we  present  the  algorithm  used  to  calculate  these  coefficients. 

Source  Coefficients 

We  first  discuss  the  equations  and  algorithms  needed  to  get  the  unit 
tetrahedron  source  coefficients.  Recalling  the  source  system  equations, 

SA  =6exp(As)Afo(Xs,Ys,Zs),  (91) 

Su=Sa[1-P(Xs,Ys,Zs)],  (92) 

Sv  =  Su -SA  ^i(Xs,Ys,Zs),  (93) 

and 

Sw  =SV -SAft2(Xs,Ys,Zs),  (94) 

where  SA ,  Sy,  Sv ,  and  Sw  are  the  cell  source  average  and  first  moments 
(known  from  the  previous  iteration),  theftj  functions  were  previously  defined 
in  Appendix  A,  and  the  coefficients  are, 

Xs=-Bu, 

Ys  =  -(Bu  +  Bv ),  (95) 

Zs  =  -(Bu  +  Bv  +  Bw). 
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This  system  of  equations  is  easily  converted  to  a  non-linear  system  of 
three  equations  and  three  unknowns  by  dividing  equations  (92)  through  (94) 

by  SA  and  letting  Puv  =-Su  ~Sy  ,  Pvw  =Sv~s w  ,  and  pw  =|w.  Thig 

fc>A  bA  bA 

substitution  and  the  application  of  Identity  6  (Appendix  A),  results  in  the 
following  system  of  equations: 

Puv  =  Pu  “  Pv  =  (Xs , Ys  ,ZS ),  (96) 

Pvw  =Pv  ~Pw  =  ^2(Xs,Ys,Zg),  (97) 

and 

pw=ft3(Xs,Ys,Zs),  (98) 

where  0  <Psum  =Puv  +  Pyw  +Pw  <1- 

The  solution  to  the  above  source  equations  is  contained  within  a  phase 
space  bounded  by  the  tetrahedron  depicted  below. 
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Figure  35.  The  (Puv>Pvw>Pw)  Solution  Space. 

The  above  phase  space  maps  (puv-pvwPw)  to  the  (XS,YS,ZS) coefficient 
space  according  to  the  mapping  below. 


Pw  0 


Figure  36.  Mapping  of  (puv,Pvw>Pw)  Phase  Space  to  (XS,YS,ZS)  Coefficient 

Space. 
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From  Figure  36,  we  see  that  the  solution  space  maps  the  input  data  in  a 
tetrahedron  within  a  volume  bound  between  0  and  1  to  all  of  (X. Yc  ,ZC ) 

space.  Additionally,  the  solutions  Xs(pUy,pyW,pw),  Ys(puy,pyW,p^),  and 
Zs(Puv  >Pvw  >Pw)>  are  expected  to  vary  rapidly  with  the  input  data  but  there 
should  be  no  discontinuities  that  would  cause  bad  behavior  with  the 
evaluation  of  the  Jacobian  that  will  be  necessary  for  the  root  solver. 

After  obtaining  Xs,  Ys,  and  Zs,  the  source  coefficients  are  found  to  be: 

Bu=-Xs, 

BV=XS-YS,  (99) 

BW=YS-ZS, 


and 


As  =loge 


Sa  \ 

v6Mo(Xs,Ys,Zs)  J 


(100) 


The  root  solve  strategy  used  in  TETRAN  to  obtain  the  source 
coefficients  is  presented  in  the  next  section. 


Source  System  Root  Solver 

The  source  coefficients  presented  above  are  obtained  by  root-solving 
the  system  of  equations  embodied  in  (96),  (97),  and  (98).  Broyden's  method 
(Burden,  1993),  which  requires  the  evaluation  of  the  Jacobian  of  the  system  of 
equations,  is  used.  We  chose  to  use  Broyden’s  method  because  it  requires  only 
one  evaluation  of  the  system’s  Jacobian  matrix  (which  is  computationally 
expensive)  followed  by  iterative  updates  that  are  computationally 
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inexpensive.  This  approach  is  taken  with  penalties  of  only  superlinear 
convergence  (versus  quadratic  for  Newton’s  method)  and  the  need  for  an 
accurate  first  guess  generator  (like  Newton's  method)  to  start  the  solver.  A 
first  guess  is  necessary.  Naively  using  (0,0,0)  as  the  first  guess,  Broyden’s 
method  failed  to  solve  the  source  system  after  20+  iterations.  Thus,  a  good 
first  guess  for  the  method  is  necessary.  However,  it  is  shown  below  that  the 
cost  to  evaluate  the  Jacobian  dominates  the  cost  of  the  root-solving  algorithm 
including  the  first  guess  algorithm. 

The  source  system  of  equations,  recast  in  the  normal  form  for  root¬ 
solving,  is 

fi(Xs,Ys,Zs)  =  ^?1(Xs,Ys,Zs)-pUy  =0,  (101) 

f2(Xs,Ys,Zs)  =  K2(Xs,Ys,Zs)-Pvw  =0,  (102) 

and  f3(Xs,Ys,Zs)  =  7?3(Xs,Ys,Zs)-Pw  =0.  (103) 


The  Jacobian  for  the  above  system  is 


rJn 

d12 

T  \ 

d13 

j  Source  _ 

J21 

J22 

d23 

d32 

d33y 

(104) 


where  the  elements  are 
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dfi 


Jll=^  =  ^(X-Ys’Zs)^l(X3’YS>Zs)-2^l(Xs,XslYs,Zs)], 

Ji2=^-=«2(X„Y„Z,)[R1(X„Y„Z.)-«1(X„Y.,Y.,Z,)], 
Jl3=^-  =  «3(X„Y„ZS)[«1(X5,YS.ZS)-«?1(XS,Y„ZS,ZS)], 

j  -  df2  -  T 

j21“ax7‘Jl2' 

J22=-^-  =  «2(Xs,Ys,Zs)[«2(Xs,Y8,Zs)-2*2(Xe,Ys,Ys,Z8)], 
J23=^--^3(XS,YS,ZS)[^2(XS)YS)ZS)-^2(XS,YS,ZS,ZS)]) 

J  -  dfs  -  J 

31“ax7“Jl3, 

J32  =  =  J23  > 


5YS 

af, 


J33=^-  =  ^3(Xs,Ys,Zs)[^3(Xs.Ys,Zs)-2^3(Xs.Ys,Zs,Zs)]. 


(105) 


In  evaluating  the  elements  of  the  Jacobian,  we  used  Identity  3  (Appendix  A) 
as  well  as  the  quotient  rule  for  differentiation. 

Clearly,  the  computation  of  the  source  system's  Jacobian  is  expensive. 
It  is  efficient  to  re-use  the  double  repeated  argument  functions  which  are 
calculated  in  the  evaluation  of  the  system  (f1;  f2 ,  and  f3)  before  the 

evaluation  of  J.  However,  this  still  leaves  6  five-argument/three-argument 
ratios  (7?jwith  repeated  arguments)  to  compute  accurately  and  robustly.  If  we 

use  the  standard  Newton's  method  to  invert  this  system,  J  would  have  to  be 
evaluated  for  each  iteration  of  the  root  solver.  Broyden's  method  (Burden, 
1993)  avoids  this  requirement  by  requiring  only  one  evaluation  of  J  followed 
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by  comparatively  cheap  iterative  updates.  Broyden  requires  1  vector  addition, 
3  matrix  multiplies,  1  transposition,  1  dot  product,  and  the  evaluation  of  fx , 

f2 ,  and  f3 .  This  approach  provides  coefficients  that  are  accurate  to  9  digits 
with  no  more  than  5  iterations  and  generally  2  or  3.  Key  to  using  Broyden's 
method  was  the  development  of  a  first  guess  algorithm,  presented  next.  The 
algorithm  used  to  numerically  evaluate  the  source  system  and  its  Jacobian  is 
presented  in  Appendix  C. 

Source  First  Guess  Algorithm 

Broyden's  Method  may  converge  slowly  or  not  at  all  if  it  is  started  too 
far  from  the  true  root  of  the  system  it  is  solving.  To  ensure  that  this  not 
happen  in  TETRAN,  we  developed  a  fairly  accurate  first  guess  algorithm 
based  upon  the  behavior  of  the  solution  space.  Indeed,  Brennan  developed  a 
similar  algorithm  based  upon  the  asymptotic  behaviors  of  Xgfpjjy  >Pw )  > 

Ys(Puv  ’Pvw  >Pw )  >  and  Zs  (Puv  >  Pvw  >  Pw )  f°r  the  solution  in  the  corners  of  the 
(Puv>Pvw>Pw)  phase  space  (Brennan,  1996).  A  simple  interpolation  spliced 
these  solutions  together,  providing  a  guess  to  the  solution  in  the  interior  of 
the  space.  This  approach  worked  well  but  did  not  take  advantage  of  the 
nature  of  the  solution  space.  Our  approach,  although  perhaps  not  optimal,  is 
based  upon  the  behavior  of  Xs (Puv  ,Pvw  ,pw ) ,  Ys(Puv  ,Pvw  ,Pw ) ,  and 

Zs(Puv  >Pvw  >Pw )  in  the  interior  of  the  solution  space. 
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Our  first  guess  strategy  is  built  around  the  asymptotic  behavior  of  the 
solution  to  the  system.  These  solutions  are  presented  in  Table  19. 

Table  19.  Asymptotic  Solutions  to  the  Source  System. 


Puv  1 

Pvw  0 

Pw  0 

PSum  1 

Puv  0 

Pvw  1 

Pw  ->  0 

PSum  1 

Puv  0 
Pvw  0 

Pw  ->1 

PSum  1 

Puv  0 

Pvw  0 

Pw  ->  0 

PSum  “*■  0 

Xs 

1 

1  1 

1  1 

1 

1  ”  PSum 

PUV  1  -  PSum 

Puv  1  -  PSun 

Puv 

Ys 

1  1 

1 

1  1 

1 

Pvw  1  -  PSum 

1  ~  PSum 

Pvw  1  ~  PSun 

Pvw 

2s 

1  1 

1  1 

1 

1 

PW  1  “  PSum 

PW  i-PSum 

1_PSum 

Pw 

In  the  case  where  any  of  Puv ,  pyw »  or  Pw  are  greater  than  or  equal  to  0.95 
and  pSum  >  0.99  or  pSum  <  0.01 ,  the  solutions  in  Table  19  are  accurate  to  at 
least  9  digits  requiring  no  root  solving. 

The  more  interesting  case  occurs  when  (pm-  ,pw  ,pw )  is  not  located  in 
one  of  the  corners  of  the  phase  space  tetrahedron.  In  this  case,  we  developed 
an  algorithm  based  upon  the  geometric  form  of  the  solutions  in  the  space. 
Using  AVS/Express,  the  solution  to  the  system  for  many  random  values  of 
(Puv  >Pvw  ’Pw )  was  visualized.  Looking  at  the  form  of  the  Xs(puy  ,py\y  ,pw ) 

solution,  for  example,  an  interesting  result  emerged.  It  is  plotted  in  Figure  37 
below. 
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Figure  37.  Plot  of  Xs(puy  ,py^,p^)  for  Xs — 100,  -10,  0,  10,  100  (right  to 

left). 

Figure  37  shows  that  Xs(puy  ,pvw  ,pw )  lies  on  constant  surfaces  that  sweep 
the  puy  axis.  These  surfaces  are  fairly  planar  for  negative  Xs  and  tend  to 
bend  (like  the  pages  of  a  book)  for  positive  Xs .  What  is  important  to  see  is 
that  Xs(puv  ,Pvw  >Pw )  appears  to  vary  in  a  one-dimensional  fashion 


as~  Xs(puv) .  Due  to  the  symmetry  of  the  system  of  equations, 

Ys(Puv-Pvw>Pw)  and  Zs(puv,pvw>Pw)  behave  in  a  similar  manner. 

The  algorithm  used  to  generate  our  first  guess  is  based  upon  the 
behavior  in  Figure  37.  Treating  the  problem  as  quasi-one-dimensional,  we 
treat  each  input  coordinate  (pw ,  pw ,  pw )  separately.  The  algorithm  is 

developed  for  the  pyy  case  because  there  is  exchange  symmetry  among  puy , 
Pvw  >  and  pw .  The  algorithm  is  as  follows. 

First,  point  data  was  generated  for  three  constant  Xs  surfaces: 

Xs  =  -5 ,  Xs  =  0 ,  and  Xs  =  5 .  These  surfaces  are  chosen  since  it  is  over  this 
range  that  the  asymptotic  behavior  breaks  down.  Each  set  of  points  was  fit 
using  TableCurve®  3D  to  generate  an  estimate  of  the  surface  equations  for 
the  three  cases.  This  was  done  once  and  put  into  a  Fortran  module.  These 
equations  and  their  statistics  are  presented  below: 

Pu iv(Pvw>Pw)  =  0.191232898  +  Pvw*  (-0.00879861  +  -0.17756583*  pyw)  + 
pw  * (-0.00919211  +  -0.1772146* pw)  +  -0.31946783* Pvw  * pw 

n  (106) 

R2  =  0.9986184864, 

Standard  Error  =  0.0021141815, 
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Puv°(Pvw>Pw)  =  0.5 -0.5*  Pvw  -0.5* pw, 

R2=1.0,  (107) 

Standard  Error  =  1.16187e- 10 

and 

Pw5(Pvw>Pw)  =  0.808632414  +  Pvw*  (-0.98752806  +0.17397106*  Pvw  )  + 

Pw* (-0.98678519  +  0.173236422* Pw)  +  0.309475303* Pm,  *Pw 
R 2  =0.9999311819, 

Standard  Error  =  0. 0020301284. 

The  above  equations  are  used  to  locate  the  input  pyy  coordinate  with  respect 
to  each  surface.  The  exchange  symmetry  is  such  that  pyw  — >  pw  and 
Pw  Puv  f°r  Ys ,  and  pw  puy  and  pw  ->  pw  for  Zs  where  puy ,  pw , 
and  pw  are  the  variables  in  equations  (106)  through  (108). 

Locating  puy  with  respect  to  the  above  surfaces,  the  value  of  Xsis 
approximated  in  one  of  four  ways.  For  pyy  to  the  left  of  the  p^)5  surface 

(Puv  —  Puv  >  0)»  Xs  = - - - .  Note  that  this  is  the  same  form  as  the 

PUV  1  “  PSum 

asymptotic  form  in  Table  19.  The  reason  for  this  is  that  the  asymptotic 
solution  is  a  very  good  approximation  in  this  range  (accurate  to  three  or  more 
digits).  When  Puyis  between  the  Xs  =  5.0  and  Xs  =  0.0  surfaces,  the 
following  fit  is  used  to  obtain  Xs : 
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xfit( Pvw)  ~  ~^°ge(P Vw) 
y  fit  (  P  Sum  )  =  ~l°§e  —  P  Sum  ) 

Xs(Pvw  ’PSum)  =1-305400076  +  Xfit*  (-4.39382576  +  Xfit*  (7.230625174  + 
xfit*  (-2.82155197  +  0.532512991*  xfit)))  + 
yfit  *(1.123135173  +  yfit  *(-4.28162193  +  (109) 

yfit  *(1.623314683  +  yflt  *(-0.31139215  + 

-0.01494551*  yfit)))), 

R2  =  0.9994254009, 

Standard  Error  =  0. 0363020055. 

In  the  case  where  pyy  is  between  the  Xs  =  0.0  and  Xs  =  -5.0  surfaces,  the  fit 
below  is  used  to  calculate  Xs : 


xfit  ( PVW  )  ~  ~l°ge  ( PVW  ) 
y fit( PSum )  ~  ~  PSu/n ) 

Xs(pvw  ,pSum)  =- 1.10819185  +  xfit*  (-1.18202335  + 

xm*  (4.431595918  +  xnt*  (-1.74308591  + 

Xfit*  (0.341900634  +  0.012610457  *  xfit))))  +  (110) 

yfit  * (3.9891 71814  +  yfit  * (-6.98404183  + 
yfit  *(2.779136284  + -0.53138969* yfit)))  , 

R2  =  0.9994254009, 

Standard  Error  =  0.0363020055. 

Finally,  if  the  value  of  pyy  is  to  the  right  of  the  Xs  =  -5.0  surface 

(Puv  -  Puv  <  0 )  the  appropriate  asymptotic  solution  (Table  19)  is  used 
because  these  solutions  are  generally  good  to  many  digits  in  this  range. 

The  source  system  root-solving  algorithm  is  summarized  as  follows. 
First,  we  evaluate  the  source  ratios  pyy ,  pw ,  and  pw  using  source  data 
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from  the  previous  iteration.  If  this  data  falls  in  the  range  of  accuracy  of  the 
asymptotic  limit  for  the  solution  to  the  system,  we  use  the  asymptotic 
solution  and  exit  the  algorithm.  Otherwise,  we  use  the  input  data  to  generate 
a  first  guess  for  Xs ,  Ys  ,  and  Zs  using  the  first  guess  algorithm  which  uses 
either  the  asymptotic  solutions  or  a  fit  to  the  data  within  several  regions  in 
the  solution  space.  Since  the  system  is  symmetric,  we  use  the  same  fit 
equations  for  each  coefficient  and  exchange  variables  in  a  consistent  way. 
After  obtaining  the  guess  for  the  coefficients,  Broyden's  method  is  used  to 
invert  the  system.  With  the  coefficients  in  hand,  we  back  out  the  needed 
coefficients  (Band  As)  and  proceed  with  the  transport  algorithm  for  the  cell. 
As  discussed  in  the  Chapter  4,  a  problem  was  found  with  the  break  point 
between  the  asymptotic  solution  and  the  root-solver.  The  solution  was  to 
change  the  break  point  to  force  root-solving  to  occur  more  often.  The  previous 
algorithms  embody  this  strategy. 

It  is  probably  possible  do  a  lot  better  than  the  above  algorithm.  The  fit 
equations,  although  picked  from  the  best  available,  could  be  better.  The 
logarithmic  versions  were  used  since  they  seem  to  better  replicate  the 
behavior  of  the  solution  of  this  exponential  system.  Additionally,  we  could 
have  broken  the  space  up  into  more  regions  and  fitted  these  with  better 
functions.  Finally,  with  more  time,  a  family  of  parametric  curves  might  have 
been  found  that  governs  the  solution  space  behavior,  perhaps  yielding  good 
enough  solutions  to  eliminate  root-solving  entirely.  Regardless,  the  above 
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guess  generator  is  good  enough  to  allow  Broyden's  method  to  converge  on  the 
average  of  2  or  3  iterations  and  no  more  than  5,  satisfying  our  need  for  a  good 
guess  generator. 

Inflow  Flux  Coefficients 

The  coefficients  for  the  inflow  flux  moments  are  obtained  in  a  manner 
similar  to  the  source  coefficients.  The  equations  for  the  inflow  flux  moments 
are 


TA=2exp(Af)AVXf,Yf), 

dll) 

=1/A[l-P(Xf,Yt)], 

(112) 

and 

'f,v='J,u-^A«i(Xf,Y(). 

(113) 

As  before  with  the  source  moments,  equations  (112)  and 

u/in  __u/in  ii/in 

(111)  and  define  p(jy  =  —  . — —  and  py  =  -4- .  Then, 

ym  q/m 

(Appendix  A)  the  face  system  becomes 


(113)  are  divided  by 
using  Identity  6 


Puv  -  Ki<Xt,Y() 

(114) 

and 

p^=«2(Xt,Yf). 

(115) 

The  above  system  is  the  two-dimensional  analog  of  the  three- 
dimensional  system  for  the  source.  It  is  non-linear  and  requires  root  solving. 

The  solution  approach  for  the  coefficients  is  the  same  as  the  source  system, 
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i.e.  use  the  asymptotic  behavior  of  Xf(p{jv,py)  and  Yf(puv,py)  as  the  basis 
for  a  first  guess  algorithm  and  employ  Broyden’s  method  as  the  root-solver. 
After  obtaining  Xf  and  Yf ,  the  needed  coefficients  are  calculated  just 

as  was  done  for  the  source  coefficients  (except  there  may  be  up  to  3  inflow 
faces).  The  coefficients  are 


B'  =-xf, 

B'=X(-Yf, 


(116) 


and 


Af  =  log€ 


Ta 


2Mo(Xf)Yf). 


(117) 


The  following  sections  present  the  inflow  flux  root-solver  equations  and 
the  first  guess  algorithm. 


Inflow  Flux  System  Root  Solver 

As  with  the  source,  we  chose  to  use  Broyden’s  method  to  invert  the 
inflow  flux  system  of  equations.  Thus,  these  equations  must  be  cast  into  the 
appropriate  format  along  with  the  system's  Jacobian. 

The  appropriately  cast  system  for  the  above  strategy  is  given  as 

f1(Xf,Yf)  =  ^1(Xf,Yf)-puv=0, 

f2(Xf  ,Yf)  =  7?2(Xf  ,Yf)-pv  =  0.  ( 


The  Jacobian  for  the  above  system,  Jf,  is 
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(119) 


Jf 


(  jf 
J11 

v4i 


j 

j 


f  ^ 
12 
f 

22  J 


where 


J 


f 

li 


jf 

d12 

f 

d21 

jf 

d22 


£ L 

axf 


^1(Xf,Yf)[^1(Xf,Yf)-2^1(Xf,Xf,Yf)]) 


£ L 

5Yf 

df2 

5Xf 

df2 

SYf 


=  (X  f ,  Yf  )[^?i  (X  f ,  Yf )  -  ^  (X  f ,  Yf ,  Yf )], 

__  jf 
"  d12> 


=  «2<Xf,Yf)[«2(Xf,Yf)-2K2(X{,Y£,Yt)]. 


(120) 


As  with  the  source  system,  the  evaluation  of  Jf  is  computationally  expensive. 
In  this  case,  it  is  efficient  to  re-use  the  three-argument  to  two-argument 
moment  function  ratios  since  they  are  needed  as  part  of  the  evaluation  of  fx 
and  f2 .  The  evaluation  of  the  repeated  argument  functions  is 

unavoidable.  Thus,  Broyden’s  method  was  chosen  to  avoid  multiple 
evaluations  of  Jf.  The  first  guess  algorithm  for  this  system  is  presented  below. 
The  algorithm  used  to  numerically  evaluate  the  inflow  face  flux  system  and 
its  Jacobian  is  presented  in  Appendix  C. 

Inflow  Flux  First  Guess  Algorithm 

The  first  guess  algorithm  for  the  inflow  face  flux  system  of  equations  is 
the  two-dimensional  analog  of  the  three-dimensional  version  used  for  the 
source  system.  Indeed,  the  inflow  face  flux  exhibits  the  same  type  of  solution 
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layering  shown  for  the  three-dimensional  version.  In  the  case  of  the  face  flux, 
the  phase  space  is  a  unit  triangle  instead  of  a  tetrahedron  (Figure  38). 


(0,1) 


Figure  38.  The  (pfjv.py)  Phase  Space. 

As  before  with  the  source  system,  the  solutions  to  the  face  system  lie 
on  well-defined  curves  of  constant  Xf(p{rv,p^)  (and  Yf(pfJV,p^)  by 
symmetry).  This  behavior  is  depicted  in  Figure  39. 
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Figure  39.  Plot  of  Xf(puy,py)  for  Xf  =-100,  -10,  0,  10,  100  (right  to  left). 

Just  as  was  observed  for  the  source,  the  above  solution  space  maps  the 
p  coordinates,  each  of  which  is  constrained  to  lie  between  0  and  1,  to  the 
entire  set  of  real  numbers  for  Xf  and  Yf . 

The  first  guess  strategy  is  essentially  the  same  as  for  the  source 
system.  The  inflow  flux  parameters  pyy  and  py  are  treated  in  a  quasi-one- 
dimensional  way  by  treating  each  coordinate  independently.  This  can  be  done 
because  the  solutions,  Xf(pfjy,py)  and  Yf (pyy  ,py) ,  seem  to  vary  more 
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strongly  with  one  coordinate  than  the  other.  The  algorithm  for  obtaining 
Xf  (Puv  ,py)  is  presented  because  the  same  formulas  apply  for  Yf  (p(jy  ,py) 

r  r 

by  letting  pyy  =  p*  in  the  equations  below. 

The  algorithm  begins  with  the  asymptotic  solutions  to  the  face  system 
in  the  cases  where  both  pfjy  and  p^  are  approaching  0  or  one  or  the  other  is 
approaching  1.  As  before,  0  <  p|um  =  p£jy  +  py  <1. 


Table  20.  Asymptotic  Solutions  for  Inflow  Face  System. 


Puv  ->■  1 

Pv  ->0 

PSum  1 

Puv  0 

Pv  -*1 

PSum  1 

Puv  0 

py  ->0 

PSum  0 

Xf 

1 

1  1 

1 

1  “  PSum 

f  f 

PUV  1  ~  PSum 

Puv 

Yf 

1  1 

1 

1 

PV  1  “  PSum 

l“PSum 

Pv 

The  aymptotic  relations  in  the  first  two  columns  of  Table  20  are  accurate  to 

»  p 

at  least  10  9  when  pgum  >  0.99.  Column  three  is  at  least  this  accurate  when 

PSum  -  0  01 .  Thus,  the  relations  in  Table  20  are  used  as  the  solutions  to  the 

face  sytem  when  they  are  applicable.  Otherwise,  the  first  guess  algorithm  is 
used. 

The  first  guess  algorithm  uses  fits  to  a  family  of  constant  Xf  curves  to 
determine  the  location  of  the  p{jy  input  coordinate  (exchange  symmetry  is 

j? 

used  to  get  pv  using  the  same  fits).  In  the  case  of  the  face  system,  we  fit  the 
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Xf  =  10,  1,  -1,  and  -5  curves  using  TableCurve®  2D.  This  was  done  once  and 
the  results  placed  in  a  Fortran  module  in  TETRAN.  These  curves  are  used  to 
bracket  pyy  in  the  solution  space.  The  equations  of  the  fit  curves  are  shown 
below: 


Pot=1°(Pv)  =  0.100241456  +  pfv  *(  -0.01458791  + 
p fv  *(0.139338267  +  pfv*(  -0.74015265  + 
p fv  *(1.114373468  +  -0.59888171* pfv )))),  (121) 

R2  =0.9999588005, 

Standard  Error  =  0.0001883807, 


Pw=1  ( Pfv )  =  0. 4185662  +  pfv*(  -0. 35932296  +  - 0. 057943 79* pfv), 

R2  =0.9999873181,  (122) 

Standard  Error  =  0.00047861 78, 


Puv=_1(P fv)  =  0.581439345  +  pfv  *(  -0.64121428  +  0.058551 053 *pfv), 

R2  =0.9999943081,  (123) 

Standard  Error  =  0. 00045081 65 , 


and 

Puv=_5(Pv)  =  0.807850555  +  pfv*  (-0.98086771  +  0.165769953*  pfv), 

R2  =0.9999370731,  (124) 

Standard  Error  =  0.002139227 . 

After  using  the  above  to  determine  where  pyy  lies  in  the  solution 
space,  we  use  a  fit  equation  appropriate  to  the  region  Puyis  in.  For  Puyto 

r  y _ ia  r  y _ i  a  r 

the  left  of  Pyy-  ,  Puy"  -  Puy  <  0 ,  we  use  the  asymptotic  solution 
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1 


1 


Xf  =  — - - — ,  which  is  accurate  to  several  digits  in  this  range.  If  pfjy 

Puv  1  -  PSum 


lies  between  10  and  pjj^.  1 ,  the  following  fit  function  is  used: 


1  1 

x  =  —f - j — , 

Puv  1  “  PSum 
1  1 

y=—s - j — , 

Pv  ^  -  PSum 

zl  =  x*(l. 2083991 38  -  0.01536180* x), 

z2 -0.5* (1.0  +  erf((y  +1.31 619805) / (42  *  2. 080453277))),  (125) 

Xf  =  -0.56584199  +  zl  +  1.024171722*  z2, 

R2  =  0.9986085942, 

Standard  Error  =  0.0995880886, 


where  erf(x)  is  the  error  function,  which  is  evaluated  efficiently  using  a 
standard  algorithm  or  a  built-in  function.  The  above  function  closely 

resembles  a  CDF  (cumulative  distribution  function)  to  which  Xf(puy  ,py)  is 
similar  as  evidenced  by  Figure  39.  If  p(jy  lies  between  p[j^=1  and  Puy=_1,  the 
fit  function  below  is  used  to  estimate  Xf : 
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X  = 


1 


1 


f  f  ? 

Pjjv  l_PSum 

1  1 

y-—f - 7 — , 

Pv  ^  —  PSum 
zl  =  0.5  *  (1.0  + 

e rf((x  +  0.5219783834650574)  7(42*5. 
z2  =  0.5  *  (1.0  + 

e  rf(fy  + 1. 34129221 4090 703)  7  (42  *  2. 9361 1 702967799))), 
Xf  =-7.313430618887315  +  13.57575241380947  *  zl  - 

3.4969015342974*  z2  +  6.489645999251979*  zl*  z2, 

R2  =  0.999932859649556, 

Standard  Error  =  0.004964030286314623, 


368613963271195))), 


(126) 


where  we  have  again  employed  the  error  function,  erf(x).  In  the  case  where 
Puy  lies  between  p[^=_1  and  Puy=~5 ,  the  last  fit  function  is  used,  which  is 
given  as  : 

1  1 

%  f  f  ’ 

P UV  1  -  PSum 

1  1 

y  =  —f - f — , 

Pv  1  ~  PSum 

zl=  x*  (1.528744066931101  +  0.06995285682140669*  x), 

z2  —  0.5*  (1.0  +  (12 

e  rf(-(y  +  4. 094407663731 77)  7(42  *  2.970 727912877624))), 

Xf  =-0.02622091422930943  +  zl  +  0.8776311335141779*  z2, 

R2  =  0.99422295436887, 

Standard  Error  =  0.09313013094057546. 
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Lastly,  when  pfjy  lies  to  the  right  of  pjj^  5 ,  the  asymptotic  behavior  is  used 

as  a  guess  to  the  solution  of  the  system:  Xf  =  -j- - ^ — . 

PUV  1-PSum 

The  inflow  face  root-solving  algorithm  is  summarized  as  follows.  First, 
the  inflow  face  flux  moment  ratios  are  computed  for  the  current  face:  p(jy 

and  py .  If  this  data  falls  within  the  region  of  applicability  of  the  asymptotic 
solutions  to  the  system,  Xf  and  Yf  are  calculated  directly  and  the  next  input 
face  is  evaluated.  Otherwise,  this  data  is  passed  to  a  first  guess  algorithm  to 
estimate  the  needed  Xf  and  Yf  using  either  the  fits  or  the  asymptotic  values. 
These  coefficients  are  used  to  start  a  Broyden's  method  root-solver.  After 
successfully  obtaining  accurate  Xf  and  Yf ,  the  needed  coefficients  are 
calculated  and  the  next  inflow  face  is  evaluated. 

The  above  approach  is  perhaps  non-optimal.  It  requires  the  evaluation 
of  erf(x) ,  which  is  efficiently  evaluated  using  standard  algorithms.  Using 
erf(x)  probably  requires  fewer  floating-point  operations  than  the  evaluation  of 
Jf  for  each  face  and  solver  iteration.  However,  this  algorithm  is  efficient  and 
ensures  the  convergence  of  the  Broyden  method  in  less  than  5  iterations  and 
generally  between  2  and  3  for  a  solution  accurate  to  at  least  9  digits.  Further 
work  could  be  done  in  this  area  to  better  take  advantage  of  the  nature  of  the 
solution  space  as  time  warrants. 
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Appendix  C:  Stable  EC  Quadrature  Formulations 

The  development  of  stable  algorithms  to  evaluate  the  flux  moments  in 
Chapter  2  was  essential  to  the  implementation  of  the  EC  method.  This 


appendix  outlines  the  algorithms  used  to  calculate  i^Jce11,  H'A.m*11,  V®usbrcce11, 


stibecll  subcell  subcell  subcell  subcell  outface 
Tu,m  7  Y  c,src  >  Yu, in  »  Y  u;,src  >  Yu;, in  >  YA,src  > 


lff  outface 

VA.in  > 


outface 
r  u,src  9 


VuUkace-  V °uSrcaCe >  an<i  ¥°Uinace-  Similar  techniques  are  used  to  evaluate  the 


source  root  solve  system  and  its  Jacobian  elements  and  the  inflow  face  flux 
system  and  its  Jacobian  elements  found  in  Appendix  B.  Note  that  throughout 
the  discussion  below  we  refer  to  the  exponential  moment  function  identities 
found  in  Appendix  A. 


IEEE®  754  Overflow  and  Underflow  Errors 

Overflow  and  underflow  errors  are  produced  because  computers  use 
finite  precision  arithmetic.  An  IEEE  overflow  is  one  where  the  resultant 

number  is  greater  than  1. 7976931348623158  x  10308  (normalized)  for  64-bit 
(double  precision)  arithmetic.  This  corresponds  approximately  to  exp (709. 78) . 

An  underflow  occurs  if  the  resultant  number  is  less  than 
2.2250738585072013  x  10~308  (normalized)  which  corresponds  to 
approximately  exp(-708.4) .  Either  case  is  intolerable  to  TETRAN.  In  the 

overflow  case,  an  error  is  produced  and  the  code  terminates.  More  insidious, 
the  underflow  case  will  gracefully  truncate  the  value  to  0.0  (or  some  other 
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small  number)  without  user  intervention.  The  EC  quadratures  require  that 
Va  >  Vw  >  V|/v  >^10  so  that  if  an  underflow  occurs  in  one  of 
the  moment  evaluations,  care  must  be  taken  to  ensure  that  the  other 
moments  correspondingly  underflow  to  0.0  as  well.  In  some  cases,  a  number 
may  underflow  to  the  un-normalized  value  of 

4.94065645841246544  x  10~324  and  not  to  0.0.  This  must  be  prevented  from 
occurring  during  the  transport  calculations. 

Source  Contribution  to  Cell  Flux  Moments 

The  source  contribution  to  the  cell  flux  moments  are  the  most 
complicated  formulae;  their  treatment  will  be  discussed  first.  In  this  case, 

'Ku^snf^ >  Vv^src11*  an<f  are  functions  of  ^j(x,y ,z , w) ,  requiring  the 

evaluation  of  multiple  five-argument  exponential  moment  functions. 
Additionally,  care  must  be  taken  to  prevent  potential  overflow  errors  because 
in  some  cases  (particularly  for  poorly  shaped  tetrahedra)  the  coefficients  will 
produce  exponents  much  greater  than  709  in  the  moment  function  routines. 
However,  it  is  known  from  empirical  evidence  that  0  <  ^j(x,y  ,z ,  w)  <  1 . 

Indeed,  it  is  the  ratio  function  7?j(x,y,z  ,w)  which  enables  the  stable  and 

robust  evaluation  of  the  EC  quadrature.  The  pseudo-code  algorithm  for  the 
evaluation  of  the  source  contribution  to  the  flux  moments  is  presented  below. 

Sort  arguments  x,  y,  z,  and  w  into  ascending  order:  the  new  list  is  xl, 
x2 ,  x3 ,  and  x4 ,  where  xl<  x2  <  x3  <  x4 .  Note  that  the  original  argument 
list  is  retained. 
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If  xl>70  0  and  as  <  -70  0  then  exp(as)is  very  small  as  are 

^l(x,y /Z  /w),^2(x/Y/z  /w)^3(x/y/z  ,w),and </?4(x/y/z  ,w)  . 

subcell  subcell  llf  subcell  lir  subcell 
VA,src  '  V u,src  '  V  v,src  •  M^w,src  0 

Else,  if  xl>  3  8.0  and  |as|<  700 then  use  the  asymptotic  value  of 
^?j(x,y,Z,w)  and  the  exp(as)will  not  overflow  or  underflow. 

ftl(x,y,z,w)  =  i,  <R2 (x,y,z,w)  =  — ,  7?3(x/y/z#w)  =  i, 
x  y  z 

ft4(x,y,z,w)  =  — 
w 


Afo(x,y,z,w)  =  ■ 


xyzw 


r^-i  subcell  subcell  llfsubcell  lirsubcell 

Calculate  ^A,src  '  Vu,src  '  Vv,src  •  M^w,src 


Else,  if  xl>  38.0  and  as>700then  an  overflow  can  occur  in  evaluating 
subcell  (  \ 

M^src  as  exP(asj  overflow.  In  this  case: 


^1(x,y,z,w)  =  l/  «2(x,y,z,w)  =  -,  K3(x,y,z,w)  =  -, 
x  y  z 


7?4(x,y,z,w)  =  —  , 
w 


A10(xfy/z,w)  = - 

xyzw 

^A.src11  =exp(as +loge(6/Mo(x,y,z,w 


Calculate 


subcell 
t  u,src  j 


subcell 
t  v,src 


...  subcell 
Vw,src  usin9 


quadrature  formulas . 


Else,  if  xl  >  -700  and  x4  <  700  and  xl-x4  >  -700  and  as  >  -700  then  use 
Identities  2  and  4  (Appendix  A)  to  evaluate  the  ^functions  using 

p functions  that  are  calculated  using  our  four-argument  moment  function 
routines . 


^1(x,y,z,w)  =  p(-x,y-x,z-x,w-x) 
^2(x,y/z,w)  =  p(x-y>-y/z-y,w-y) 
ft3(x,y,z,w)  =  p(x-z,y-z  ,-z,w-z) 
^4(x,y,z,w)  =  p(x-w,y-w,z-w,-w) 
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Calculate 


subcell 
rA,src  ' 


subcell 
T  u,src  ' 


subcell  subcell 
t  v,src  *  T  w,src 


Else,  if  xl  >  -700  and  x4  <  700  and  xl-x4  >  -700  and  as  <  -700  then  use 
Identities  2  and  4  (Appendix  A)  to  evaluate  the  functions  using 
p functions  that  are  calculated  using  our  four -argument  moment  function 

routines.  The  exponential,  exp(as)  ,  will  underflow  so  use  the 
logarithmic  formulation  presented  previously  for  M^srcP  • 


^l(x,y,z,w)  =  p(-x  ,y-x,z-x,w-x) 
■R2(x,y,z,w)  =  p(x-y  ,-y,z-y,w-y) 


^3(x,y,z,w)  =  p(x-z,y-z,-z,w-z) 
<R4(x,y,z,w)  =  p(x- w,y- w,z- w,-w) 

V  A  S11  =  exP(as  +  loge(6ZMo(x,y,z,  w)» 


Calculate 


subcell  subcell 
' r  u,src  •  T  v,src 


subcell 
•  t  w,src 


using  quadrature  formulas . 


Else,  use  a  stable  IEEE  754  compliant  routine  to  evaluate 

‘Z?1(x,y/z,w),^2(x,y,z,w),^3(x,y,z,w),and^4(x,y/z,w) .  in  this 
algorithm,  we  apply  the  appropriate  set  of  identities  for  the  input 

argument  order  of  which  there  are  24  different  cases.  This  is  because 

the  argument  order  is  sorted  in  ascending  order  but  the  original 
arguments  must  be  known  in  order  to  use  the  correct  algorithm.  Thus, 
the  sorted  argument  lists  could  be  one  of  the  following:  xyzw,  xywz, 
xzyw,  xzwy,  xwyz,  xwzy,  yxzw,  yxwz,  yzxw,  yzwx,  ywxz,  ywzx,  zxyw,  zxwy, 

zyxw,  zywx,  zwxy,  zwyx,  wxyz,  wxzy,  wyxz ,  wyzx,  wzxy,  and  wzyx.  It  is 

from  one  of  the  above  argument  lists  that  the  needed  functions  are 
calculated.  The  example  algorithm  for  the  xyzw  case  is  presented  below. 
This  algorithm  can  be  permuted  for  any  of  the  above  given  cases. 

Using  Identities  2  and  4  (Appendix  A) : 


!Z?1(xfy,zfw)  =  p(-x,y-x,z-x,w-x) 

Using  Identities  1  and  2  (Appendix  A) : 

M0(-x/y-x,y-x,z-x)-A1o(y- 


^2(x,y,z,w)  = 


Alo(-x,y-x,z-x)-  M)(y- 


x,y  -  x,z  -  x,w  -  x) 
x,z  -  x,w  -  x) 


7?3(xfy,z,w) 


A\q(-x,y  -  x,z  -  x,z-  x)-  A%(y-  x,z  -  x,z  -  x,w  -  x) 
At0(-x,y-x>z-x)-Aio(y-x,z-x,w-x) 


7?4(x,y,z,w) 


A1o(-x,y  -  x,z  -  x,w  -  x)-  Mq{y  -  x,z-x,w-x,w-x) 
A1o(-x,y-  x,z-  x)- A(o(y-x,z-  x,w  -  x) 


Using  Identity  2  (Appendix  A) : 
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. , ,  subcell 
M;A,src 


6/exp(as  -x)A1o(~x,y  -  x,z  -  x,w  -  x) 


Calculate 


subcell 
t  u,src  ' 


llf  subcell  subcell 

V  v,src  -  V  w,src  using  quadrature 


formulas . 


The  above  IEEE  754  compliant  algorithm  uses  Identity  1  to  avoid  the  need 


for  a  general  five-argument  exponential  moment  function  routine  (  Mq  and 


p).  If  such  a  routine  were  available,  the  above  algorithm  would  be  much  less 


complicated  because  one  would  only  worry  about  which  coeffiecient  was  xl. 


Inflow  Flux  Contribution  to  Cell  Flux  Moments 

The  algorithm  for  calculating  v|/^eU,  v|/®u?nceU,  x\ /®u?nceU,  and  \|/®u^eU  is 

analogous  to  the  algorithm  above  except  that  7?j(x,y,z)  is  the  needed 
function. 

Sort  arguments  x,  y,and  z  into  ascending  order:  the  new  list  is  xl,  x2 , 
and  x3 ,  where  xl<  x2  <  x3  .  Note  that  the  original  argument  list  is 
retained. 

If  xl>700  and  af  <  -700  then  exp(af)is  very  small  as  are 

<^l(x,y,z),^2(x,y/z),and^3(x,y,z)  . 

...  subcell  .ir  subcell  subcell  %l, subcell 

VA.in  '  Vu,in  f  Vv,in  *  Vw,in  are  0 

Else,  if  xl>  38.0  and  |af|<700then  use  the  asymptotic  value  of 
7?j(x,y ,z)  and  the  exp(af)will  not  overflow  or  underflow. 

<R1(x,y,z)  =  —  l  <R2(x,y,z)  =  -,  ^3(x,y,z)  =  l 
x  y  z 


Mo(x,y,z)  = 


xyz 


r>=i  .-...i  ...subcell  subcell  subcell  .,,subcell 
Calculate  l|/A>in  ,  V|/uin  ,  V|/v>in  ,  ywin 
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Else,  if  xl>  3  8.0  and  Rf  >700  then  care  must  be  taken  in  evaluating 
subcell  /  \ 

M^A,in  as  exp^afj  will  overflow.  In  this  case,  use  an  alternate  but 
equivalent  form. 


^1(x,yfz)  =  l,  ft2(x, y,z)  =  ~,  'R3(x,y,z)  =  - 

X  y  Z 


Ato(x,y,z)  = - 

xyz 


VAUineU  =  exp(af  +loge(6Ato(x,y,z))) 


Calculate 


subcell 
•  u,in  ' 


subcell 
tv, in  • 


subcell 
t  w,in 


using  quadrature  formulas. 


Else,  if  xl  >  -700  and  x3  <  700  and  xl-x3  >  -700  and  af  >  -700  then  use 

Identities  2  and  4  (Appendix  A)  to  evaluate  the  functions  using 

p functions  that  are  calculated  using  the  three -argument  moment  function 
routines . 


^!(x,y,z)  =  p(-x,y-x,z-x) 
7?2(x,y,z)  =  p(x-y,-y,z-y) 
ft3(x,y,z)  =  p(x-z,y-z;-z) 


Calculate 


.(r  subcell 

VA.in  - 


subcell 
t  u,in  ' 


subcell  subcell 
t  v,in  >  t  w,in 


Else,  if  xl  >  -700  and  x3  <  700  and  xl-x3  >  -700  and  af  <  -700  then  use 
Identities  2  and  4  (Appendix  A)  to  evaluate  the  ^functions  using 
p functions  that  are  calculated  using  three -argument  moment  function 
routines.  The  exponential,  exp(af)  ,  will  underflow  so  use  the 

Sll  nPP!  1 

logarithmic  formulation  presented  previously  for  \\f  jn 

^l(x,y,z)  =  p(— x,y-x,z-x) 


7?2(x,y,z)  =  p(x-y,-y,z-y) 
7?3(x,y,z)  =  p(x-z,y-z  ,-z) 

V  AUineU  =  ®xP(a  f  +  loge  ( 6  Mq  (x ,  y ,  z))) 


Calculate 


subcell 
Tu,in  ’ 


subcell 
t  v,in 


subcell 
t  w,in 


using  quadrature  formulas. 


Else,  use  a  stable  IEEE  754  compliant  routine  to  evaluate 

!R1(xfy,z),^2(xfy,z),and^3(xfy,z)  .  In  this  algorithm,  use  the 
appropriate  set  of  identities  for  the  input  argument  order  of  which  we 
are  only  interested  in  three  cases:  xl  =  x,  xl  =  y,  or  xl  =  z.  This  is 
different  from  above  because  one  had  to  know  the  argument  order  to  use 
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Identity  1  effectively.  Here,  because  we  have  a  general  f our-argument 
moment  function  routine,  no  concerned  is  required  about  the  order 
following  xl .  The  example  algorithm  for  xl  =  x  is  present  below.  This 
algorithm  can  be  permuted  for  the  other  two  cases. 

Using  Identities  2  and  4  (Appendix  A) : 

7?l(x,y#z)  =  p(-xfy-x,z-x) 

Using  Identities  1  and  2  (Appendix  A) : 


^2(x,y,z) 


yVto(-x,y-x,y-x,z-x) 

Mo(-x,y-x,z-x) 


^3(x,y,z)  = 


A^(-x,y-x,z-x,z-x) 


Afo(-x,y-x,z-x) 

Using  Identity  2  (Appendix  A) : 

^subceU  =6eXp(af  -x)Afo(-x,y-x,z-x) 


Calculate 


subcell 
i  u,in  • 


subcell 
T  v,in  • 


subcell 
t  w,in 


using  quadrature  formulas. 


Source  Contribution  to  Outflow  Face  Flux  Moments 


The  algorithm  used  to  determine  v|/^Ug^ce, 


outface 
t  u,src 


,  and  \\i 


outface  • 
t»,src 


identical  to  that  for  the  inflow  flux  contribution  to  the  cell  flux  moments.  The 


only  difference  is  the  requirement  for  three  flux  moments  and  the  equation 

forvSS*. 


Sort  arguments  x,  y,and  z  into  ascending  order:  the  new  list  is  xl,  x2 , 
and  x3 ,  where  xl<  x2  <  x3  .  Note  that  the  original  argument  list  is 
retained. 


If  xl>700  and  as  <  -700  then  exp(as)is  very  small  as  are 
^l(xfy,z),^2(x,y,z),and^3(x,yfz) . 


lt  outface 
VA.src 


outface 
t  u,src 


and 


outface 
t  v,src 


are  0 


Else,  if  xl>  38.0  and  |asj<700then  use  the  asymptotic  value  of 
^j(x,y,z)  and  the  exp(as)will  not  overflow  or  underflow. 


<Rl(x,y,z)  =  —  ,  (R2(x,y,z)  =  -l  <R3(x,y,z)  =  - 
x  y  z 
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Afo(x,y,z)  = 


1 


xyz 


Calculate 


outface  outface 
rA,src  '  Tu?src  • 


and 


outface 
T  v,src 


Else,  if  xl>  38.0  and  ag  >700 then  care  must  be  taken  in  evaluating 
outface  /  \ 

V  A,src  as  exPlasj  wil1  overflow.  In  this  case,  use  an  alternate  but 
equivalent  form. 


ftl(x,yfz)  =  i,  ^2(x,y,z)  =  l>  7?3(x,y,z )  =  - 
xyz 


Aio(x,y,z)  = 


,  outface 
VA.src 


xyz 

exp(as  +loge(2/Mo(x,y,z))) 


Calculate  '  an<^  V|/^^ce using  quadrature  formulas. 

Else,  if  xl  >  -700  and  x3  <  700  and  xl-x3  >  -700  and  as  >  -700  then  use 
Identities  2  and  4  (Appendix  A)  to  evaluate  the  ^functions  using 

p functions  that  are  calculated  using  the  three -argument  moment  function 
routines . 


(x,y,z)  =  p(-x,y-x,z-x) 


ft2(x,y,z)  =  p(x-y,-y,z-y) 
7?3(x,y,z)  =  p(x-z,y-z  ,-z) 


Calculate 


outface 
V  A,src  - 


outface 
t  u,src 


and 


outface 
t  v,src 


Else,  if  xl  >  -700  and  x3  <  700  and  xl-x3  >  -700  and  as  <  -700  then  use 
Identities  2  and  4  (Appendix  A)  to  evaluate  the  ^functions  using 
p functions  that  are  calculated  using  the  three -argument  moment  function 

routines.  The  exponential,  exp(as) ,  will  underflow  so  use  the 

outface 

logarithmic  formulation  presented  previously  for  M^A^rc  * 


<R1(x.,y,z)  =  p(-x,y-x,z-x) 


7?2(x,y;z)  =  p(x  — y  ,-y,z-y) 


7?3(xfy,z)  =  p(x-z,y-z,-z) 
vSSrc58  =  exp(as  +loge(2Z  Mo(x,y,z))) 


Calculate 


outface 
t  u,src 


and  \)/°Ug^ce using  quadrature  formulas. 
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Else,  use  a  stable  IEEE  754  compliant  routine  to  evaluate 

^l(x/y/z),^2(x/y/z)^nd^3(x/y,z)  .  In  this  algorithm,  apply  the 
appropriate  set  of  identities  for  the  input  argument  order  of  which  we 
are  interested  in  only  three  cases:  xl  =  x,  xl  =  y,  or  xl  =  z .  This  is 
different  from  above  because  one  had  to  know  the  argument  order  to  use 
Identity  1  effectively.  Here,  because  there  is  a  general  four-argument 
moment  function  routine,  there  is  no  concern  about  the  order  following 
xl.  The  example  algorithm  for  xl  =  x  is  present  below.  This  algorithm 
can  be  permuted  for  the  other  two  cases. 

Using  Identities  2  and  4  (Appendix  A) : 


^1(x/y,z)  =  p(-xfy-x,z-x) 

Using  Identities  1  and  2  (Appendix  A) : 


ft2(x/ y,z)  = 


Mp(-x,y  -  x,y-x,z-x) 
M)(-x,y-x,z-x) 


Mo(-x,y -x,z-x,z-x) 

<R3(x,y,z)  =  — *--■ - - - 

Alo(-x,y-xfz-x) 

Using  Identity  2  (Appendix  A) : 

vS66  =  21  exp(as  -x)Ato(-x,y  -  x,z  -  x) 


Calculate 


outface 
■  u,src 


and  Vvs^rc06 using  quadrature  formulas. 


Inflow  Flux  Contribution  to  the  Outflow  Face  Flux  Moments 

As  before,  this  case  is  just  a  lower  dimensional  version  of  the  preceding 
cases.  This  time,  the  quadrature  formulas  depend  on  .  The  flux 


moments  calculated  are  > 


outface 
•  u,in 


,  and  \|/ 


outface 

u,in 


Sort  arguments  x,  and  y  into  ascending  order:  the  new  list  is  xl  and 
x2 ,  where  xl<  x2  .  Note  that  the  original  argument  list  is  retained. 


If  xl>700  and  af  <  -700  then  exp(af)is  very  small  as  are 
^j(x/y) and7?2(x,y)  . 


, , ,  outface 
A, in  ' 


outface 
t  u,in 


and 


outface 
t  v,in 


are 


0 


Else,  if  xl>  38.0  and  |af|<700then  we  can  use  the  asymptotic  value  of 
7?j(x,y)  and  the  exp(af)will  not  overflow  or  underflow. 
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Ki(x,y)  =  — ,  ^2(x/Y)  =  - 

x  y 


Mo(x,y)  =  — 
xy 


Calculate 


. ,  r  outface 
VA.in  ' 


outface 
t  u,in  • 


and 


outface 
t  v,in 


Else,  if  xl>  38.0  and  af  >700then  care  must  be  taken  in  evaluating 
outface  /  \ 

VA,in  as  exP(afJ  will  overflow.  In  this  case,  use  an  alternate  but 
equivalent  form. 


«1(xfy)  =  l,  Kjj(x, y)  =  - 

x  y 


Mo(x,y)  =  — 

xy 

^outfare  _  exp(af  +loge(2M0(x,y))) 


Calculate 


outface 

Yu,in 


,  outface 
and  l|/v#in  using 


quadrature  formulas . 


Else,  if  xl  >  -700  and  x2  <  700  and  xl-x2  >  -700  and  af  >  -700  then  use 
Identities  2  and  4  (Appendix  A)  to  evaluate  the  ^functions  using 

p functions  that  are  calculated  using  our  two-argument  moment  function 
routines . 


^i(x,y)  =  p(-x,y-x) 

^2(x/y)  =  p(x-y -y) 


Calculate 


...  outface 
V  A, in  ' 


outface 
T  u,in  < 


and  V \) 


outface 

v,in 


Else,  if  xl  >  -700  and  x2  <  700  and  xl-x2  >  -700  and  af  <  -700  then  use 
Identities  2  and  4  (Appendix  A)  to  evaluate  the  ^functions  using 
p functions,  which  are  calculated  using  the  two-argument  moment  function 

routines.  The  exponential,  exp(af) ,  will  underflow  so  use  the 
logarithmic  formulation  presented  previously  for  . 


^i(x/y)  =  p(~x/y-x) 


^2(x/y)  =  p(x-y/-y) 

y  outface  _  exp(af  +loge(2Mo(x,y))) 


Calculate 


outface 
T  u,in 


,  ...  outface 
and  l|/  v  jn  using 


quadrature  formulas . 
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Else,  use  a  stable  IEEE  754  compliant  routine  to  evaluate 

^?1(x/y)and<R2(x/Y)  *  In  this  algorithm,  apply  the  appropriate  set  of 
identities  for  the  input  argument  order  of  which  we  are  only  interested 
in  two  cases:  xl  =  x  or  xl  =  y.  The  example  algorithm  for  xl  =  x  is 
presented  below.  This  algorithm  can  be  permuted  for  the  other  case. 

Using  Identities  2  and  4  (Appendix  A) : 

^i(x,y)  =  p(-x>y-x) 

Using  Identities  1  and  2  (Appendix  A) : 


^2(x/Y)  = 


Mo(-x,y-x,y-x) 


wio(-x»y-x) 

Using  Identity  2  (Appendix  A) : 

^ outfaee  _  2  exp(af  _  x)A1o(-X,y  -  x) 


Calculate  and  v[/®u^ace  using  quadrature  formulas. 


Source  System  and  Jacobian  Elements 

For  the  sake  of  brevity,  we  will  not  present  the  algorithm  for 

evaluating  source  system  functions.  The  source  system  is  comprised  of 

z)  functions  whose  evaluation  was  previously  discussed  in  the  section, 

Source  Contribution  to  Outflow  Face  Flux  Moments. 

As  for  the  Jacobian  elements,  an  approach  similar  to  that  discussed  in 
Source  Contribution  to  Cell  Flux  Moments  was  used.  In  this  case,  we  have 
five-argument  moment  functions  over  three-argument  function  ratios  that 
are  functions  of  only  three  distinct  arguments  with  more  than  one  repeat 
coefficient.  These  functions  are  products  of  functions.  Briefly,  the  needed 

Jacobian  elements  are: 
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j  1 1  =  ^1  (x,  X  z)[R 1  (x,  y,z)~  2  (x,  X,  y,  z)], 

Ji2  =  ^2  (x  X  «)[^i  (x  X «)  -  ^1  (XX  y,  z)], 
J13  =  ^(XX  z)[R. i  (X  y,  z)-n 1  ( x,  y,  z,z)\ 
^21  =  Jl2> 

J22  =  ^(XXZ^CXX*)  -2 ^(XXX*)]. 

J23  = ' ^3  (x  y>  z)[^2  (x>  x z)  -  ^2  (x>y>  z> «)]. 

J3I  =  Jl3> 

J32  =  *^23  > 

J33  =  ^(XX^I^XX*)  -  2<Rz(x,y,z,z)\ 


To  evaluate  the  Jacobian  elements  that  are  products  of  functions,  we  use 

an  algorithm  similar  to  that  used  in  the  section,  Source  Contribution  to  Cell 
Flux  Moments.  The  pseudo-code  for  this  algorithm  is  shown  below. 


Sort  arguments  x,  y,and  2  into  ascending  order:  the  new  list  is  xl,  x2, 
and  x3 ,  where  xl<  x2  <  x3  .  Note  that  the  original  argument  list  is 
retained. 

If  xl  >  38.0  then  we  use  the  asymptotic  solutions  for  the  needed  ratios 
presented  below. 


^i(x,y,z)^i(x,x,y/z) 


<R2{-x.,y,-z)<Rl(yi,ylylz) 


^3(x,y,z)«i(x,y,2,z) 


7?2(x,y,z)7?2(x,y,y,z) 


!R3(x,y,z)ft2(x,y,z,z) 


<Rz{x.,y,z)1i3{yi,ylz,z) 


Alo(x,x 

,x,y, 

z) 

1 

M)(x 

.Y,z) 

M)(X/X 

'Y'Y, 

z) 

1 

Mq(x 

'Y'Z) 

xy 

A1o(x'x 

,Y,z, 

z)  _ 

1 

Mq(x, 

rY/Z) 

xz 

Mq(x,  y 

•  Y,Y, 

z) 

1 

M)(x 

.Y.z) 

=7 

M)(x»y 

,Y.z, 

z)  _ 

1 

M)(x 

,Y,z) 

yz 

M)(x/ y- 

fz  ,** 

z)  _ 

1 

Ato(x, 

y,z) 
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Else,  if  max[abs(xl),  abs(x3)]  <  1.0  then  apply  Identities  2  and  4 
(Appendix  A)  to  the  numerator  to  find  the  Jacobian  elements. 


Rl(XiYtz)Rl(x,x,ylz) 


Ato(x,x,x,y,z) 

Ato(x,y,z) 


=  exp(-x) 

Mo(x,y,z) 


4R2(X'Y'*)llk(x,Y,Y,*) 


A-to(x,x,y,y,z) 

Ato(x,y,z) 


=  exp(-x) 


M(~x>y~x/y~x>z~x) 

Mo(x,y,z) 


^3(x,y,z)!Ri(x,y,z,z) 


Mq(x, x,y,z,z) 
Mq(x, y,z) 


=  exp(-x) 


A4i(-x;y-x,z-x,z-x) 

A1o(x,y,z) 


ft2(x,y,z)7?2(x,y,y,z) 


M)(x,y,y,y,z)  =  Mi(x~y,  ~Y,  0,z-y) 

M)(x.y/Z)  Mo(x,y,z) 


^3(x/y/z)^2(x,yfz,z) 


Ato(x,y,y,Z>z)  _  ,  |M(x~y>  -y,z-y,z-y) 

A1o(x,y,z)  A<o(x,y,z) 


^3(xfy,z)^3(x,y,z,z) 


Mq(x ,y ,z  fz  ,z) 
Ato(x,y,z) 


=  exp(-z) 


Al^x-z^-z,  -z,  0) 
Ato(x,y,z) 


Else,  xl<  38.0  but  the  arguments  are  not  small.  In  this  case,  use  an 
algorithm  similar  to  the  source  contribution  algorithm  in  which  the 
sort  list  is  important.  Fortunately,  there  are  only  six  cases:  xyz, 
xzy,  yxz ,  yzx,  zxy,  and  zyx.  We  present  the  case  where  the  sort  list  is 
xyz.  The  other  cases  are  permutations  on  this  case.  Note  that 

•Ali  =  AIq  p  .  The  previous  case  was  needed  because  the  following 
algorithm  is  unstable  for  all  coefficients  small  and  close  together 
(Identity  1) . 

Applying  Identities  2  and  4  (Appendix  A) : 


<R1(xly,z)<Rl(x,x,y,7:) 


AjoC^x^x, y/z) 
Alo(x,y,z) 


M1(-x,  0,y-x,z-x) 
Mq  (-x  ,y  -  x  ,z  -  x) 


^(x.y.z^ix.y.y.z) 


Mp(x, x,y,y,z)  ^A\{-x,y-x,y-x,z-x) 
Mo(x,y,z)  Afo(-x,y-x,z-x) 


^3(x,y/z)!Ri(x/y,zfz) 


Mo(x,x,y,z,z)  _Mi(-xly-x,z-x,z-x) 
Mo(x,y,z)  Mo(-x,y-x,z-x) 


Applying  Identities  1  and  2  (Appendix  A) : 
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ft2(x,y,z)ft2(x,y,y,z) 


Mo(x,y,y,y,z) 

Mo(x,y,z) 


^H)(~x,y-x/y-x<y-x)- At0(y-x>y-x<y-x,z-x) 
M0(-xfy-x)-M0(y-x,z-x) 


ft3(x,y,z)ft2(x,y,z 


AjgQc^y^y^z) 

Xo(x,y,z) 

_  >H)(~x>y-x,y-x,z-x)- Ato(y-x,y-x,z 
Xq  (-x  ,y  -  x)  -  A1o(y  -  x  ,z  -  x) 


x,z-x) 


ft3(x,y,z)ft3(x,y,z 


_  Mo(x,y,z,z,z) 

Mo(x,y,z) 

Mo(-x,y-x,z-x,z-x)- Ato(y-x,z-x,z 
M)(-x,y  -  x)  -  Ato(y  -  x,z  -  x) 


x,z-x) 


The  above  algorithm  would  be  more  efficient  if  a  general  five-argument 
moment  function  routine  existed.  However,  the  algorithm  is  stable,  robust 
and  accurate  without  such  a  routine. 


Inflow  Face  Flux  System  and  Jacobian  Elements 

Similar  to  the  source  system,  the  inflow  flux  system  is  composed  of 
ftjO.y)  functions,  which  are  evaluated  using  the  same  approach  used  in 

Inflow  Flux  Contribution  to  the  Outflow  Face  Flux  Moments.  Thus,  we  will 
not  discuss  this  algorithm  further. 

As  with  the  source  system  Jacobian,  the  inflow  face  flux  Jacobian 
elements  are  comprised  of  products  (x,y)  and  products  of  functions. 

Briefly,  these  Jacobian  elements  are: 
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Jn  =  <Ri{x>y)[^i{x,y)-2(R1{x,x,y)\ 
Ji2  =  fy&y^i&y)-R\(x,y,yj\, 

jf  _  jf 
d21  ~  d  12  > 

J22  =  <^2(x>y)[<^2ix>y)-2'^2(x>y>y)\ 


The  algorithm  for  evaluating  j/)has  been  presented.  The  algorithm  used 
to  evaluate  the  other  ^  products  is  analagous  to  that  used  for  the  source 
Jacobian  and  is  shown  below. 


Sort  arguments  x  and  y  into  ascending  order:  the  new  list  is  xl  and  x2 
where  xl<  x2 .  Note  that  the  original  argument  list  is  retained. 

If  xl  >  38.0  then  use  the  asymptotic  solutions  for  the  needed  ratios 
presented  below. 

^(x,y)«l(x,x,y )-^;X|,‘y)-4 

M)(x^y)  x2 

^.y)*,(X,y,y)  =  ^^=i 

Mo(x/Y)  xY 

R2(x,y)«2(x,y,Y)  =  ^p^  =  4 

Ato(x,y)  y2 

Else,  if  xl<  38.0  use  Identities  2  and  4  along  with  stable  routines  to 

calculate  two-,  three-,  and  four-argument  moment  functions.  We  present 
the  case  of  xl  =  x.  The  other  case  mirrors  this  case.  Note  that 

Mi  =Mq  p  . 

Applying  Identities  2  and  4  (Appendix  A) : 

K1(x,y)K1(x,x,y)  =  °-y-f 

Mo(x,y)  Mo(-x,y-x) 

M)(x/ y)  Mo(-x,y-x) 

R2(x,y)R2(x,y,y) = ^x;y,y:y)=yH,(~x'y;x,y~x'y~x) 

Mo(x,y)  A1o(-x,y-x) 
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Appendix  D:  Derivation  of  the  Surface  Cell  Algorithm 
Motivation 

We  began  the  development  the  surface  cell  spatial  quadrature  to 
address  a  mesh  problem  that  is  quite  difficult  for  all  current  unstructured 
mesh  generators.  Imagine  a  simple  cube  surrounded  by  a  thin  wall  l/100th 
the  thickness  of  the  cube  (or  the  thin  skin  of  an  aircraft  or  satellite  or  55 
gallon  drum  of  biological  or  chemical  agent).  When  we  mesh  this  problem 
using  tetrahedra,  the  mesh  generator  will  do  one  of  three  things.  It  might 
produce  a  coarse  mesh  with  several  poor  aspect  ratio  tetrahedra  in  the  thin 
region.  Or,  it  might  attempt  to  fill  the  thin  region  with  many  (thousands)  of 
small,  well-shaped  tetrahedra.  Or,  finally,  it  might  crash  because  of  the 
memory  requirements  for  making  a  proper  mesh  for  this  type  of  problem.  The 
impact  of  the  first  case  is  that  the  poorly  shaped  tetrahedra  are  not 
numerically  well  conditioned  and  will  introduce  spatial  and  numerical  errors 
into  the  transport  problem.  The  second  case  is  difficult  in  that  the  small  mesh 
dimension  will  propagate  to  areas  of  the  mesh  where  we  can  tolerate  a  coarse 
mesh,  significantly  driving  up  the  computational  cost  of  the  problem.  The 
final  case  prevents  the  solution  of  the  problem. 

These  mesh  problems  stem  from  the  requirement  (for  transport)  that 
each  tetrahedron  in  a  mesh  share  the  faces  of  its  neighbors.  We  require  that 
the  faces  be  shared  because  the  flux  moments  from  an  upstream  cell  are 
passed  to  downstream  cells  via  affine  transformations  of  the  moments  on  the 
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shared  faces.  Presumably,  finite  element  algorithms  for  stress  and  thermal 
require  connected  meshes  for  similar  reasons.  The  requirement  for  connected 
cells  is  a  difficult  constraint  to  satisfy  with  a  minimal  number  of  well  shaped 
cells  when  the  interfacing  regions  of  a  mesh  differ  substantially  in  geometric 
extent,  i.e.  thin  regions  next  to  thick  regions. 

Tetrahedron  Aspect  Ratio 

E 

The  aspect  ratio  for  a  tetrahedron  cell  is  defined  as  y  =  — ,  where  E  is 

h 

the  length  of  the  longest  edge  and  h  is  the  shortest  height  (distance  from  a 
node  to  opposite  face).  This  ratio  is  shown  in  Figure  40. 


Figure  40.  Definition  of  Tetrahedron  Cell  Aspect  Ratio. 

We  can  see  from  Figure  40  that  if  the  cell  is  smashed  flat  or  stretched  long, 
the  aspect  ratio  will  be  large.  We  seek  to  eliminate  these  types  of  cells  by 
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using  surface  cells.  The  minimum  possible  aspect  ratio,  1.2247,  is  achieved  by 
using  an  equilateral  (regular)  tetrahedron. 

Returning  to  our  cube  example,  Figure  41  shows  the  first  case  (poor 
aspect  ratio  tetrahedra)  for  the  simple  cube  surrounded  by  a  thin  layer  of 
material.  The  mesh  on  the  left  contains  42  tetrahedra  and  the  mesh  on  the 
right  contains  4902  tetrahedra. 


Figure  41.  Examples  of  Unstructured  Tetrahedra  Meshes  with  Thin  Regions. 

In  the  coarse  mesh  case  (the  left  cube),  36  out  of  the  42  tetrahedra  (all  in  the 
shield)  have  aspect  ratios  of  144!  Figure  42  shows  a  histogram  of  the 
distribution  of  aspect  ratios  for  the  finer  mesh  on  the  right.  Note  also  that  the 
right  mesh  is  quite  irregular  around  the  corners  of  the  cube,  presumably 
attempting  to  detail  possible  stress  risers  for  a  structural  finite  element 
model. 
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0.40 


Figure  42.  Distribution  of  Cell  Aspect  Ratios  for  Fine  Meshed  Example 

Problem. 

Obviously,  the  distribution  of  aspect  ratios  is  better  than  in  the  coarse  mesh 
case  for  Figure  42.  However,  there  are  still  an  overwhelming  fraction  of  cells 
with  large  aspect  ratios  in  this  mesh.  Indeed,  the  only  way  to  dramatically 
reduce  the  average  aspect  ratio,  other  than  using  surface  cells,  is  to  mesh  the 
problem  with  cells  whose  maximum  dimension  is  on  the  order  of  the  shield 
thickness.  This  mesh  is  not  shown  because  in  every  attempt  to  produce  such 
a  mesh  I  locked  up  the  mesh  generator,  Parametric  Technology  Corporation's 
Pro/Mesh™,  Release  18.0  (Parametric  Technology  Corporation,  97).  The 
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hardware  used  in  these  attempts  was  the  ASC  MSRC's  SGI  Power  Onyx, 
which  has  1  Gigabyte  of  memory  (ASC,  1998)! 

The  Surface  Cell  Algorithm 

Clearly,  there  is  need  for  better  ways  to  approach  thin  volume  regions 
in  an  unstructured  mesh.  Our  approach  is  the  surface  cell  algorithm. 
Consider  Figure  43  below. 


Figure  43.  Cell  and  Side  View  of  a  Surface  Cell. 

Figure  43  shows  an  extruded  triangular  prism  along  the  local  eT  axis  of  the 
( eu,ev,ew )  coordinate  system.  This  cell  has  an  intrinsic  thickness,  Aw,  and 

an  optical  path  length  of  with  respect  to  the  streaming  direction,  Q . 


Particles  enter  the  front  face  at  w  =  0  and  exit  the  back  face  at  w  =  Aw . 
Actually,  a  few  particles  enter  or  exit  through  the  sides  but  we  will  ignore 
these  since  we  assume  that  Au;  is  small  compared  to  the  lateral  extents  of 
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the  cell.  We  erect  the  (eH  ,et,,eu,)  coordinate  system  to  span  the  surface  cell  as 
though  it  were  a  triagular  prism.  Indeed,  the  surface  cell  is  reprented  as 
having  no  thickness  in  thhe  mesh  (one  would  have  to  take  measures  to 
preserve  the  surface  cells  material  in  implementing  the  quadrature). 

The  key  approximation  for  the  surface  cell  algorithm  is  that  we  can 
neglect  the  lateral  displacement  of  the  flux  exiting  the  cell  and  treat  the 
transport  in  the  cell  in  a  slab-like  manner  (Figure  44).  Clearly,  as  Aw  -»  0, 
this  approximation  becomes  increasingly  accurate. 


Figure  44.  Surface  Cell  Approximation. 

Note  that  in  Figure  44  we  have  increased  the  path  length  across  the  cell  to 
match  the  path  le'ngth  of  the  correctly  propagated  ray  (dashed  arrow). 

With  the  geometric  foundation  set  and  approximations  presented,  we 
are  now  ready  to  derive  the  surface  cell  spatial  quadrature  for  the  linear  and 
exponential  characteristic  methods. 
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Derivation  of  Surface  Cell  Spatial  Quadratures 

As  with  characteristic  methods  based  on  a  tetrahedron  cell,  the  surface 
cell  method  is  also  based  on  the  characteristic  equation  for  the  cell  flux, 

\\)(u,v,w)  =  f(u,n)exp(-^)  +1^-  S(u,v,w')exp(— :SW~W  2)  (128) 

A  w  J0  ^  Aw 

where  f(u,u)  is  a  distribution  for  the  inflow  face  flux  and  %  is  the  extended 
path  length  across  the  surface  cell.  As  will  be  shown,  f  (u,v)  need  not  actually 

be  specified  ,  nor  any  form  assumed.  Only  its  moments  will  be  needed.  The 
outflow  flux  is  the  flux  at  w  =  Aw , 

Aw;  i  r  r 

\\fout(u,v)  =  \\f(u,v,Aw)=f(u,v)exp(-e)  +  exp(-e)  f— S (u,v,w')  exp(^~).  (129) 

o  £  Aw 


Lastly,  we  define  a  cell  volume  moment  operator  for  the  surface  cell, 


M[g]  =  2  J  duj  du  J  ^g(u,v,w), 

Ann 


(130) 


where  this  is  the  volume  operator  for  a  regular  triangular  prism.  If  we 
separate  the  lateral  and  axial  integration  operators,  we  can  define  M  as 

M[g(u,  v,w)]  =  L[  g(u,v)]  A[g(ic)], 

A[g(m)]  =  [-^g(w),  (131) 

0  A  w 

1  u 

L[  g(u,v)  ]  =  2  J  duj  d  vg(u,v). 

o  o 
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Indeed,  the  key  approximation  for  the  surface  cell  algorithm  is  the  separation 
of  variables  as  implied  by  the  operators  in  (131). 

We  now  derive  the  linear  and  exponential  characteristic  versions  of  the 
surface  cell  algorithm. 

Linear  Characteristic  Surface  Cell  Approximation 

We  begin  the  linear  characteristic  surface  cell  derivation  by  assuming 
that  the  normalized  source  is  distributed  linearly  in  the  la-dimension  across 
the  cell, 

S(u,  v,  w) = [P0(iu) + 6  Pj  (w)  ]  g(u,  v),  (132) 


where  P0(iu)  =  1,  0  =  and  g (u,v)  is  a  source  distribution 

Aw  SA 

whose  average  is 


1  U 


Aw  j  , 

M[ S(u, v, ie)] = 2  J du  J du g(u, v)  j  — — [P0 (u/)  +  0  Pj  (zu')  ] 


0  o 

1  U 


Aw 


(133) 


=  2jdujdv  g(u,v)=SA . 


0  0 


Thus,  the  functional  form  of  the  cell  flux  is 


\j/(u, v,w)  =  f(u, v ) exp(-s  — ) + fdu/[l+0P1(u/)]  exp[  ^].  (134) 

Aw  g  J0  Aw 
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For  the  outflow  flux,  we  start  with  (129)  and  substitute  the  source 

distribution,  noting  that  2 - 1=1 -2(1 - ) .  Changing  variables  such  that 

A  w  Aw 

wf 

t= - ,  the  outflow  flux  becomes 

Aw 


v out  (u,v)=f(u,u)  exp  (-s)  +  g(u,  V )  exp(-s)  —  K0  (-s) + 

*  (135) 
0  g(u,  V )  exp  (-8)  ■ [  2 Kx  (-s)  ^ -  K0  (-s )] , 

where  the  K  -functions  are  functions  originally  presented  in  the  linear 
characteristic  paper  (Mathews,  1998)  and  are  defined  as 


0  0  0 


(136) 


Operating  on  (134)  with  A,  we  calculate  the  zeroth  moment  of  the  flux 
with  respect  to  w, 

\]/0(u,v)=A[\\i(u,v,w')] 


(137) 


g(u,v)Acdw 


6"tdw'wr  -z(w'-w”) 

l-—|du;"[1+eP1(^")]exp( — - -). 

o Aw o  Aw 


IX) f  XX)'' 

Changing  variables  such  that  t'  = - and  t"  = - ,  we  transform  (137)  into 

Aw  Aw 
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(138) 


1 

Vo(u>v)  =  Jdr  f(u,u)exp(-st')  + 


g(u,v)Aw 


i  r 


JdrJdr[i+e(2r-i)]exp(-s(*'-r)). 

0  0 


Changing  variables  again  such  that  t  =  V  - 1"  and  performing  the 
integrations,  we  find  that 


Aw 


V  o  ( u,  v )  =  f  (u,v)K0  (e)  +  g(u,  v)—K00  (e)  + 

§ 

e  g (U.U)  ^[2Xu(s)-2XM(e)-ir0i0(8)]. 


(139) 


With\|/0(u,L>)  in  hand,  we  now  derive  the  first  moment  of  the  flux  with 
respect  to  w,  i|/j(u,v) .  As  before,  we  are  only  evaluating  the  inner  integral  of 
M  using  A.  In  the  first  moment  case, 


\\i1(u,v)=A[3P1(w')\\i(u,v,w')]. 


(140) 


Using  the  same  changes  of  variables  as  for  \\i0(u,v) ,  we  find  that 


Aw, 


^i(u)v)=3i(u>v)[2K1(s)-K0(z)]+3g(u>v)~[2Kw(s)-K00(s)]+ 


Aw, 


(141) 


3  0  g(u,v)—[4K20(z)~  4Ku(e)-  4K10(s)  +  2 K0,(e)+  K0fi(e) J. 


With  the  above  moments,  we  are  ready  to  define  the  surface  cell 
approximation.  We  define  the  results  of  the  following  moment  operator 
operations: 
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(142) 


Ui(u,v)]=y™, 

L[uf(u,u)]=M/;r, 

L[uf(u,u)]  =vj/jn, 

L[g(u,u)]=Syl, 

L[ug(u,i>)]  =  S1( , 

L[ug(u,u)]=Sl,, 

recalling  that  we  defined  the  source  normalization  to  ensure  that 
l[g(u,v)]=SA .  The  cell  flux  moments  based  on  the  source  normalization  are: 


Aw 


M/f  =L[vi/0(u,o)]=m/“^0(s)+S^—  K0'0(e)  + 


QSA^[2Kh0(s)-2K0A(e)  -Kofi(e)], 


(143) 


Aw 


v|/“u  =  L[uvj/0(u,o)]=\i/“  K0(e)+Su  —K00( s)  + 


A  w. 


0S„  ^[2Kl0(s)-2K0l(s)  -K0S)(b)], 

5 


(144) 


Aw 


iiCU  =L [v\\i0(u,v)]=vT  K0(e)  +  Sv—K00(e)  + 


0  St,— - —  [2i^10(s)  —  2iC01(s)  iC00(s)], 


(145) 


and 


Aw, 


v|/“n  =L[  v|/1(u,u)]=3\|/“[2K1(£)-X0(s)]  +  ^^/i  "Y"[2  K10(e)  Z£00(s)]  + 


3  0  SA^-[ 4  ^2,0  (e)  - 4  Ku  (e)  -  4  tf1>0 (e)  +  2  K0A  (e)  +  K0,0  (s)] . 


(146) 


171 


Operating  on  \|/out(w,i;)  with  L  gives  us  the  outflow  flux  moments: 


V T  -■ U ¥ °ut (u, y)]=v|/ “  exp(-e) +SA  ^K0( s)  + 
0S  a^(K0(s)-2K1(s)), 


(147) 


Aw 


mC  =  L[w  \\i out  (u,  u)] = \[/  “  exp(-s)+Su  —Kq(z)  + 

% 

es„  ^-(K,(z)-2K,(t)), 


(148) 


Aw 


MC  =L[av(/0Ut(u,i;)]=ii/“exp(-e)+Sy—  K0(e)  + 

0S„  ^-(K0(z)-2K,(z)). 


(149) 


Finally,  we  derive  the  conservation  relationships  by  applying  M  to  the 
BTE,  just  as  we  did  in  the  tetrahedron  cell  case.  The  conservation  equations 


are: 


Average: 


vf-vJ+svr-S, 


Aw 

T! 


(150) 


u-moment: 


V? -v* +*V?  =8, 


Aw 

T’ 


(151) 


o-moment: 


vr-v;+s<"=s^, 


(152) 


and  ia-moment:  3  \\i  +  3  \\i  “  -  6  \\)  ^e11  +  s  v|/  ®u  =  S  ^ 


w 


(153) 
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The  step  characteristic  quadrature  is  obtained  by  using  equations 
(143),  (147),  and  (150)  for  the  average  flux  in  the  cell.  With  the  linear 
characteristic  surface  cell  equations  developed,  we  now  extend  the  algorithm 
to  the  exponential  characteristic  method. 


Exponential  Characteristic  Surface  Cell  Approximation 

The  exponential  characteristic  surface  cell  assumes  an  exponentially 
distributed  source  in  the  w  direction, 


exp(-p-^) 

S (u,v,w)  = — -  g(u,v), 

Mo(P) 


(154) 


where  again  we  have  introduced  a  normalization  ( - )  such  that 

AU P) 

A[S(u;)]  =  1  and  L[g(ii,u)]  =S^ .  However,  in  this  case  we  must  determine  the 
value  of  the  source  parameter,  p,  via  root  solving  as  is  done  in  the  exponential 

characteristic  method  for  slab  geometry  (Mathews,  1993).  Indeed,  the  source 
parameter,  p ,  is  determined  by  the  root-solving  algorithm  used  for  the  EC 

method  on  unstructured  grids  of  triangles  (Mathews,  1997).  In  this  instance, 

g 

we  take  the  ratio  of  0  =  — ^  and  solve  for  p(P) ,  giving  us 

SA 

P(P)=-^,  (155) 

b 


173 


which  is  solved  for  p  using  the  method  referenced  above.  Note  that  -3  <  0  <  3 
so  that  0  <  p(P)  <1. 

We  begin  the  derivation  of  the  cell  flux  with  (128)  and  substitute  for 
S  (u,v,w) , 

\\>(u,v,w)  =  f(u,u)exp(-^)  +  f dm' exp(-p— ) exp( ~£ — ~  —  -).  (156) 

Aw  ^AIq(P)^  Aw  Aw 


For  the  outflow  flux,  we  substitute  (154)  into  (129): 


i|/out(u,  v)  =  \\i(u,v,Aw)  =f(u,i>)exp(-e)  + 


taw  exp(-p^)  exp(i^^).  <157) 

V  lAu>'  a,„ 


Aw 


As  in  the  LC  case,  if  we  change  variables  such  that  t= and  use  the 

Aw 

definition  of  the  moment  functions,  we  find  that  the  outflow  flux  is 

Vout (U, V) = f  (u, V ) exp(-s) + Aw  exp(-s)  (p  -  s).  (158) 

M)(P)  % 

As  before,  we  will  first  apply  the  A  operator  on  equations  (156)  and 
(158)  to  obtain  the  w  contribution  to  the  flux  moment.  For  the  sake  of  brevity 
and  since  the  operations  and  changes  of  variables  are  identical  to  the  LC  case 
(we  just  have  moment  functions  instead  of  K  functions),  I  will  just  present  the 
results  of  the  operations: 
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V o  ( U ,  v) = A[v(it, v,w')] = f (u, v )  (s)  +  Aw  Mp  (s,P), 

M)(P)  5 


(159) 


V 1  (m»  v) = A[  3  Pj  ( w ')  \\ f(u,  v,w ')] 

=  3f(u,if)A10(s)[l-2p(£)]  +3^^^AUs,p)[l-2p(£,»]. 

M>(P)  4 


(160) 


With  v|/0(u,u),\|/1(u,i;)  ,  and  \j/out(w, a) defined,  we  now  operate  on  them 
with  L  and  use  the  moment  definitions  in  (142),  producing  the  flux  moments 
for  the  cell: 


vf  =Uv0M] 

=v!inM.(e)+s^ 


Aw  Mq(8,P) 

4  M,( P)  ’ 


MCU  =  L[uy0(u,u)] 

=  vt/‘nM0(e)+SIt 


Aw  Mq(s,P) 

%  M>(P)  ’ 


M C“  =L[wv|/0(n,i;)] 

=  VpnM)(e)+St, 


Aw  Mp(s,p) 

5  M>(P)  ’ 


(161) 


(162) 


(163) 


and 


vf  =yVl(u,£)] 

=3y"M,(£)[l-2p(c)]+3S^A^(^>[l-2p(e,p)]. 

%  M>(P) 


(164) 


Similarly,  the  outflow  face  flux  moments  are: 
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=urut(u,v)] 

=  M/!i  exp(-e)  +S^  exp(-s)— 

%  M,(P) 

V"  =L[wvf/out(u,i;)] 

=  MC  exp(-e)  +S  exp(-s)— 

4  M)(P) 


(165) 


(166) 


and 


¥ 


out 

V 


=  L[v^t(u,v)] 

=  ¥“  exp(-e)  +SV  exp(-e)—  ^o^~s\ 

§  Mo(p) 


(167) 


The  conservation  equations  for  the  average-,  u-,  and  i>-moments  are 
exactly  the  same  as  those  for  the  LC  case;  equations  (150),  (151),  and  (152). 
The  equation  for  m-moment  conservation  is 


3mC  +  3  -6V?  +eWf=3SA  ^[l-2p(P)].  (168) 

Clearly,  we  can  see  that  if  we  can  trade  a  surface  cell  calculation  for  a 
tetrahedron  calculation,  we  gain  great  efficiencies.  The  surface  cell  algorithm 
requires  no  cell  splitting  and  in  the  case  of  EC,  only  a  single  one-dimensional 
root  solve  [which  was  optimized  by  Mathews  (Mathews,  1997)].  Both  EC  and 
LC  require  rotation  of  i|/“  and  vj/“  into  the  surface  cell  coordinate  system, 
and  determination  of  the  orientation  of  the  streaming  direction  with  respect 
to  the  cell  zero  coordinate  to  correctly  account  for  the  w  moments  as  is 
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normally  done  for  slab  codes  (by  sign  reversal).  Contrast  this  with  a 
tetrahedron  cell,  which  requires: 

-  cell  splitting, 

-  translation/rotation  of  inflow  and  source  moments, 

-  multi-dimensional  non-linear  root  solves  for  the  source  and  inflow  flux, 

-  many  calls  to  exponential  functions  for  the  spatial  quadrature, 

-  and,  re-assembling  of  the  sub-cell  moments. 

Further,  we  avoid  poorly  shaped  cells,  which  propagate  throughout  the  mesh. 

Implementation  of  the  Surface  Cell  Quadrature 

Unfortunately,  the  above  spatial  quadrature  was  not  implemented  in 
this  research.  Before  implementing  the  surface  cell  algorithm,  we  need  a 
mesh  generator  that  adequately  addresses  the  special  needs  of  a  mixed  mesh 
transport  code:  connected  meshes  that  conserve  volumes  and  materials. 
Pro/Mesh  was  inadequate  to  this  task.  Thus,  we  derived  the  quadrature  for 
future  use  and  proceeded  to  develop  the  parallel  version  of  TETRAN. 
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